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ABSTRACT 

We present a comprehensive model of the global properties of Alfven waves in the solar 
atmosphere and the fast solar wind. Linear non-WKB wave transport equations are solved from 
the photosphere to a distance past the orbit of the Earth, and for wave periods ranging from 3 
seconds to 3 days. We derive a radially varying power spectrum of kinetic and magnetic energy 
fluctuations for waves propagating in both directions along a superradially expanding magnetic 
flux tube. This work differs from previous models in three major ways. (1) In the chromosphere 
and low corona, the successive merging of flux tubes on granular and supergranular scales is 
described using a two-dimensional magnetostatic model of a network element. Below a critical 
flux-tube merging height the waves are modeled as thin-tube kink modes, and we assume that 
all of the kink-mode wave energy is transformed into volume-filling Alfven waves above the 
merging height. (2) The frequency power spectrum of horizontal motions is specified only at 
the photosphere, based on prior analyses of G-band bright point kinematics. Everywhere else in 
the model the amplitudes of outward and inward propagating waves are computed with no free 
parameters. We find that the wave amplitudes in the corona agree well with off-limb nonthermal 
line-width constraints. (3) Nonlinear turbulent damping is applied to the results of the linear 
model using a phenomenological energy loss term. A single choice for the normalization of the 
turbulent outer-scale length produces both the right amount of damping at large distances (to 
agree with in situ measurements) and the right amount of heating in the extended corona (to 
agree with empirically constrained solar wind acceleration models). In the corona, the modeled 
heating rate differs by more than an order of magnitude from a rate based on isotropic Kolmogorov 
turbulence. 

Subject headings: MHD — solar wind — Sun: atmospheric motions — Sun: corona — turbulence 
— waves 



1. Introduction 

Magnetic fields are known to play a significant role in determining the equilibrium state of the plasma 
in the solar atmosphere and solar wind (e.g., Parker 1975, 1991; Narain and Ulmschneider, 1990, 1996; 
Priest 1999). Much of the magnetic flux in the "quiet" photosphere seems to be concentrated into small 
(100-200 km) intergranular flux tubes. The physical processes that heat the chromosphere and corona 
have not yet been identifled deflnitively, but there is little doubt that magnetohydrodynamic (MHD) effects 
are prevalent. Even many proposed nonmagnetic mechanisms depend on the underlying properties of the 
magnetically structured atmosphere. The outflowing solar wind is fed by open magnetic flux tubes, and 
many MHD processes have been proposed to deposit heat and momentum at locations ranging from the 
extended corona to interplanetary space. 
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The continually evolving convection below the photosphere; gives rise to a wide spectrum of MHD 
fluctuations in the magnetic atmosphere and wind. The propagation of waves through the solar atmosphere 
has been studied for more than a half century (Alfven 1947; Schwarzschild 1948; Biermann 1948). Although 
compressible (e.g., acoustic and magnetoacoustic) MHD waves are likely to be dynamically and energetically 
important in some regions of the atmosphere, it is the mainly incompressible Alfven mode that has been 
believed for many years to be dominant in regions that are open to the heliosphere (e.g., Osterbrock 1961; 
Kuperus, lonson, & Spicer 1981). Indeed, the MHD fluctuations measured by spacecraft in the solar wind 
have a strongly Alfvenic character (Belcher & Davis 1971; Hollweg 1975; Tu & Marsch 1995; Goldstein, 
Roberts, & Matthaeus 1995). 

Even though much has bec;n learned about the generation, propagation, reflection, and damping of 
Alfven waves in the solar atmosphere, most earlier studies have focused on only a finite range of heights and 
treated the interactions with regions above or below as boundary conditions. This necessarily involved the 
approximation that the relevant phenomena are mainly local, i.e., that they do not depend on the conditions 
very far away from the region being modeled. There are circumstances, however, for which this approximation 
breaks down. For example, the properties of reflecting Alfven waves with long periods (i.e., of order 1 day) in 
the solar wind depend formally on the conditions "at inflnity," since they behave asymptotically as standing 
waves (e.g., Heinemann & Olbert 1980). 

In this paper we construct a comprehensive model of the radially evolving properties of Alfvenic fluctua- 
tions in a representative open magnetic region of the solar atmosphere and fast solar wind. The model takes 
account of nonlocal interactions by tracing the waves from their origin as transverse flux-tube oscillations in 
the photosphere all the way out to the interplanetary medium (truncated for convenience at 4 AU). This is 
done with the smallest possible number of free parameters. There are two overall aims of this paper: 

1. We wish to understand better the global "energy budget" of Alfven waves, including relative amplitudes 
of inward and outward propagating waves, along the open flux tubes that feed the solar wind. 

2. In order to determine how MHD turbulence contributes to the heating of the extended solar corona, 
we need to pin down precisely how Alfven waves provide the natural preconditions for a turbulent 
cascade. 

The second aim above was motivated by a recent study of the small-scale dissipation of MHD turbulence in 
the extended corona (Cranmer & van Ballegooijen 2003). This kinetic dissipation depends strongly on how 
the turbulence is excited at its largest scales, and in this paper we attempt to put the "outer scale" wave 
dynamics on firmer footing so that the resulting "inner scale" can be better understood. The work described 
by this paper builds on prior studies by Hollweg (1973, 1978a, 1981, 1990), Heinemann & Olbert (1980), 
Spruit (1981, 1984), An et al. (1990), Barkhudarov (1991), Velli (1993), Lou (1993, 1994), Lou & Rosner 
(1994), MacGregor & Charbonneau (1994), Kudoh & Shibata (1999), Matthaeus et al. (1999), Hasan et al. 
(2003), and many others. 

One unique aspect of this paper is that the photospheric spectrum of transverse fluctuations — which 
constrains the Alfven wave amplitudes at all larger radii — is specifled directly from detailed observations of 
magnetic bright point (MBP) motions and is not (as in many other models) given as an arbitrary bound- 
ary amplitude. However, a complete physical description of the fluctuations in the photosphere (e.g., the 
relative inward/outward amplitudes and the kinetic/magnetic energy partition) is obtained only after the 
fully nonlocal wave reflection has been computed for all radii. Another way this work differs from many 
previous models is that the expansion and successive merging of flux tubes on granular and supergranular 
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scales is described using a two-dimensional model of a magnetic network element in the inhomogeneous solar 

atmosphere. 

Despite the attempted comprehensiveness of this model, we needed to make three specific approximations 
in order to render the calculations tractable. These approximations are summarized here, but they are also 
discussed further below and justified for certain regimes of applicability. First, we ignore all effects of 
compressible fluctuations (e.g., acoustic waves; fast and slow magnetosonic waves) despite their importance 
in understanding observed intensity oscillations and chromospheric heating. This is done in order to fairly 
assess the relative importance of the incompressible Alfven mode before resorting to more involved models. 
Second, the Alfven wave model is mainly linear, which limits its applicability in regions where the wave 
amplitudes become large in comparison to background field strengths and characteristic speeds. (Some 
effects of nonlinearity are examined, though, in § 6.) Third, we do not model explicitly the back-reaction 
of the waves on the mean properties of the solar atmosphere and solar wind. We do, however, compute 
quantities such as the wave pressure acceleration and turbulent heating rate for use in future models of this 
back-reaction. 

The remainder of this paper is organized as follows, fn § 2 we give an overview of the physical processes 
to be incorporated in the wave models together with a "cartoon" description of the steady-state magnetic 
field topology. In § 3 we describe the adopted steady-state (i.e., zero-order) plasma conditions in detail. The 
specification of the photospheric frequency spectrum of transverse fluctuations is given in § 4, and the wave 
equations to be solved arc given in § 5. Solutions of these equations, including some with various prescriptions 
for damping, are presented in § 6. We conclude with a summary of major results (§ 7) and a discussion 
of remaining issues (§8). Appendix A contains supplementary equations describing analytic solutions of 
the kink-mode wave equations for isolated flux tubes in an isothermal atmosphere. Appendix B compares 
various published formalisms for the non-WKB transport equations for Alfven waves in an accelerating 
wind. Appendix C summarizes the properties of the fully developed anisotropic MHD turbulence spectrum 
discussed by, e.g., Cranmer & van Ballegooijen (2003). 



2. Overall Picture of Open Field Regions 

Figure 1 illustrates the basic magnetic field geometry that we believe is representative of flux tubes 
that feed the high-speed solar wind. A key feature of the adopted configuration is the successive merging of 
strong-field magnetic flux tubes between granules (Figure la) and supergranules (Figure lb). On the largest 

scales. Figure Ic shows the more or less axisymmctric magnetic field that is characteristic of solar minimum 
(using the field model of Banaszkiewicz, Axford, & McKcnzic 1998), but nearly all of the work presented in 
this paper can also be applied straightforwardly to open-field regions at other phases of the solar cycle. 

We assert that most of the plasma that eventually becomes the time-steady solar wind originates in 
intergranular magnetic flux tubes known observationally as G-band bright points, network bright points, or 
in groups as "solar filigree." This assertion is seemingly uncontroversial (knowing what we do about solar 
magnetic fields), though it is seldom stated. Adopting a ripening convention in nomenclature, we refer to 
these 100-200 km size photospheric features as magnetic bright points (MBPs).^ High-resolution observations 
reveal the presence of MBPs in the dark lanes between granules, and these features are associated with regions 



^We also note that MBPs arc not the same phenomena as the larger bright points in the chromosphere (Rutten & 

Uitenbroek 1991) or the still larger X-ray bright points in the low corona (e.g., Golub et al. 1977; Parnell 2002). 
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of strong (1 2 kG) magnetic field believed to be contained within nearly vertical flux tubes (e.g., Sheeley 
1967; Dunn & Zirker 1973; Muller 1983, 1985; Piddington 1978; Rabin 1992; Solanki 1993; Berger & Title 
2001). There is increasing evidence for magnetic structures on even smaller scales than 100-200 km, but 
we leave the study of these structures to future work (see, e.g.. Stein & Nordlund 2002; Sanchez Almeida, 
Emonet, & Cattaneo 2003; Trujillo Bueno, Shchukina, & Asensio Ramos 2004). 

As magnetic flux rises stochastically from the convection zone to the photosphere (e.g.. Priest, Hey- 
vaerts, & Title 2002), field lines are simultaneously jostled horizontally by fluid motions on granular (1-2 
Mm) scales. Magnetic flux is concentrated into thin tubes by some combination of "flux expulsion" from 
the upflowing granule centers to the downflowing lanes (Parker 1963), the rapid evacuation of these down- 
flowing superadiabatic regions (i.e., convective collapse; Parker 1978; Spruit 1979), and enhanced radiative 
cooling leading to thermal relaxation (Sanchez Almeida 2001). MBPs are observed frequently to merge with 
neighboring flux elements and to spontaneously fragment into several pieces (e.g., Berger et al. 1998). Once 
formed, MBPs continue to be shaken back and forth by the underlying convective motions (Kulsrud 1955; 
Osterbrock 1961; van Ballegooijen 1986; Huang, Musielak, & Ulmschncider 1995), which results in various 
kinds of wavelike fluctuations describable in the "thin-tube" MHD limit (Spruit 1981, 1982, 1984; see further 
references in § 5.1 below). Oscillatory motions can also be induced by impulsive reconnection events (e.g., 
Moore et al. 1991) or, conversely, the random wave trains may display observational time signatures that 
could be misinterpreted as small-scale flaring (Moriyasu et al. 2004). 

There is a great variety of MHD wave activity expected and observed in the inhomogeneous solar 
atmosphere. In addition to isolated MBP fluctuations, there is ambient acoustic wave energy excited by con- 
vection, shock steepening in the chromosphere and transition region, and both driven and free oscillations in 
sunspots and coronal loops (see recent reviews by Axford et al. 1999; Roberts 2000; Bogdan 2000; Hirzberger 
2003). A signiflcant fraction of the Sun's magnetic flux may also be distributed outside the MBPs, network, 
and active regions (e.g., Schrijver & Title 2003), and acoustic waves in "field-free" regions may in fact be 
magnetoacoustic. Damping of MHD waves and turbulence has been a key ingredient in many proposed 
models of chromospheric and coronal heating. Our main focus, though, is on the incompressible waves that 
eventually escape from the atmosphere into the solar wind. 

Magnetic flux tubes rooted in MBPs undergo both transverse (kink-mode) and longitudinal (sausage- 
mode) oscillations that can propagate upward from the photosphere. Because the strong-field flux tubes 
are in horizontal pressure equilibrium with the surrounding weak-field material, they have a lower density 
and thus are susceptible to buoyancy effects and evanescence for long enough periods. Like acoustic waves, 
the compressible longitudinal modes steepen into shocks and damp over a few scale heights, while the 
incompressible kink modes can propagate into the corona relatively undamped. Nonlinear effects can lead 
to mode conversion between kink and longitudinal modes (e.g., Ulmschneider, Zahringer, & Musielak 1991). 

Somewhere in the low chromosphere, the thin flux tubes are believed to expand laterally to the point 

where they merge with one another into a homogeneous network field distribution (Spruit 1984; Pneuman, 
Solanki, & Stenfio 1986; Gu et al. 1997). At this point, the thin-tube description breaks down and standard 
MHD wave theory becomes more applicable (i.e., kink modes become transverse Alfven waves). The merged 
network flux bundles have horizontal scale lengths of 2 to 6 Mm and are probably maintained by large-scale 
convective flows that push the fleld to the edges of supergranular cells. At a larger height — still below the 
chromosphere-corona transition region — the network magnetic field expands laterally to 10-30 Mm scales 
and is thought to merge again into a large-scale "canopy" (e.g., Kopp & Kuperus 1968; Gabriel 1976; 
Giovanelli 1980; Anzer & Galloway 1983; Dowdy, Rabin, & Moore 1986). The spatial scale of the canopy 
is set by the typical distance between network flux bundles (i.e., the size of supergranulation cells) in the 
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chromosphere. Observational evidence for preferential wind acceleration in the rapidly expanding network 
"funnels" is growing (Rottman, Orrall, & KHmchuk 1982; Hassler et al. 1999; Peter & Judge 1999; Aiouaz et 
al. 2004; see also Martmez-Galarce et al. 2003) but is still not definitive (e.g., Dupree, Penn, & Jones 1996). 

As waves propagate upward into the corona, the radially varying Alfven speed allows for gradual linear 
reflection (Ferraro & Plumpton 1958). The transition region can also act a sharp "reflection barrier" to 
Alfven waves with wavelengths exceeding the local scale length of the Alfven speed in that thin zone (see 
§6). It is thus possible for a time-steady STiperposition of upward and downward propagating Alfven waves 
to be maintained (e.g., HoUweg 1981, 1984). Strong reflection is not necessarily an impediment to there 
being a substantial upward wave flux; one merely needs more power in the upward modes than in the 
downward modes. Somewhere in the solar atmosphere the MHD fluctuations become turbulent, but it is 
unclear whether the turbulent cascade becomes energetically important in the photosphere (Petrovay 2001), 
in the chromosphere and transition region (Chae, Schiihle, & Lemaire 1998), or in the corona (e.g., Dmitruk 
& Matthaeus 2003). 

In the extended corona, the high-speed solar wind begins to accelerate supersonically above a heliocentric 
distance of 2 to 3 solar radii (Rq). In the radially inhomogeneous wind, the dissipationless propagation of 
MHD waves does work on the mean fluid and provides an added wave-pressure acceleration (e.g.. Belcher 
1971; Jacques 1977; Leer, Holzer, & Fla 1982). The presence of the wind also modifies how waves propagate, 
and above the Alfven critical point — where the wind speed equals the local Alfven speed at about 10 Rq — 
both the inward and outward modes are advected outward with the wind. Large coronal holes are the most 
probable source regions for the fast solar wind (for a review of observations and theoretical models, see 
Cranmer 2002). Some flux tubes in coronal holes have a higher density than the surrounding open-field 
plasma; these "polar plumes" seem to trace out the superradial expansion of the merged-canopy magnetic 
field in the corona (e.g., DeForest, Lamy, & Llebaria 2001). 

In addition to Alfven waves, there is some evidence for both fast-mode and slow-mode magnetosonic 
waves in the corona (Ofman, Nakariakov, & DeForest 1999; Nakariakov et al. 2004), but they have been 
observed only in relatively confined regions such as loops and plumes. Fast and slow modes are believed to be 
more attenuated by coUisional damping processes than Alfven waves before they reach the corona. However, 
once some fraction of the energy flux of Alfven (and possibly fast-mode) waves escapes into the solar wind, 
the classical transport theory of coUisional damping begins to break down and coUisionless wave-particle 
interactions should dominate the damping (e.g., Hollweg & Isenberg 2002; Marsch, Axford, & McKenzie 
2003; Cranmer 2000, 2001, 2002, 2004; and many references therein). Studies of wave-particle resonance 
damping have seen a recent resurgence of interest because of their potential importance in producing the 
preferential ion heating and acceleration seen both in the extended corona (Kohl et al. 1997, 1998, 1999) 
and in situ (Marsch 1999). 

As Alfven waves propagate into interplanetary space, their velocity and magnetic-field amplitudes grow 
to nonlinear magnitudes (i.e., 6B/Bq becomes of order unity). Higher order ponderomotive effects begin to 
dominate the wave propagation equations (Lau & Siregar 1996) and coUisionless wave-wave resonances can 
create additional mode coupling and damping (Lee & Volk 1973; Lacombe & Mangeney 1980; Lou 1993). At 
some point it may even be inappropriate to use a "wavelike" paradigm to discuss the increasingly turbulent 
fluctuations (e.g., Goldstein et al. 1995). In the outer heliospherc, the description of magnetic flux tubes 
as relatively closed systems — implicit in the above discussion — breaks down as well, since processes such as 
stream-stream interactions, interstellar pickup ion injection, and cosmic ray transport increasingly dominate 
the physics. 
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In this paper we purposefully examine only a subset of the many processes listed above. This is done 
mainly to keep the modeling tractable, but it also helps clarify the extent to which the specific wave modes 
that we study can account for various observations. 

3. Steady-State Plasma Conditions 

Models of linear waves depend sensitively on the assumed background zero-order plasma properties (e.g., 
density, flow speed, magnetic field strength, flux tube geometry). In this section we describe the empirically 
constrained time-steady plasma along the central axis of a radially pointed (but superradially expanding) 
magnetic flux tube, from the photosphere to a distance of 4 AU. The empirical description consists of two 
parts: a two-dimensional magnetostatic model of thin MBP flux tubes that expand into the supergranular 
network canopy (§3.1) and a one-dimensional analytic continuation of the plasma parameters in the extended 
corona, assumed here to be along the axis of symmetry of a polar coronal hole at solar minimum (§ 3.2). 

3.1. Network Magnetic Structure 

We flrst develop a two-dimensional numerical model of a supergranular network element as a collection 
of thin flux tubes. The gas pressure in the atmosphere decreases with increasing height, causing a lateral 
expansion of the flux tubes. Neighboring MBP flux tubes within the network element merge into a monolithic 
structure at some height Zm- (Heights z are measured from the optical-depth- unity photosphere; radii r are 
measured from Sun center.) For a thin, isothermal flux tube in pressure equilibrium with its surroundings, 
the interior magnetic field strength varies with height as e~^/^^, where H is the pressure scale height (e.g., 
Spruit 1981, 1984). Thus we determine the so-called merging height by solving for 

Zra « 2H\ii{Bq/B) , (1) 

where Bq is the field strength inside the flux tube at the base of the photosphere (z = 0) and B is the 
average flux density in the network patch. Observations of the MBP flux-tube field strength range between 
1000 and 2000 G, and the network-averaged field strength varies between about 20 and 300 G (e.g., Gabriel 
1976; GiovancUi 1980). In order to determine a representative value for the merging height we use H — 120 
km. Bq = 1500 G, and i? = 100 G to obtain Zm « 600 km, a height in the low chromosphere (see also 
Pneuman et al. 1986). Above the merging height the network element consists of a single thick flux tube 
that further expands with height. The outer edge of this tube forms a magnetic canopy that overlies the 
neighboring supergranular cells. A second merging occurs when neighboring network elements come together 
at a "canopy height" Zc above the supergranular cell centers; we set Zc = I Mm (see also Hasan et al. 2003). 
Figure 2 shows the three-part structure of the network element: 

1. The region below the merging height (0 < z < 0.6 Mm) is described as a collection of thin flux tubes 
embedded in a fleld-free medium. The field strength B{z) is assumed to be the same for all flux tubes, 
and their cross sections and other plasma properties are consistent with the thin-tube approximation. 

2. Between heights of 0.6 and 1 Mm, the merged network fiux element expands laterally to the edges of 
the supergranular cell, which overlies a field-free cell-center chromosphere. 

3. Between 1 and 12 Mm, the "fully merged" magnetic field fills the supergranular cell volume and expands 
primarily in the vertical direction. 
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In the remainder of this paper, the term "merging height" refers specifically to the merging of the thin flux 

tubes at Zm = 0.6 Mm. 

The total magnetic flux $ of the network element is constrained empirically to be 3.7 x 10^^ Mx. The 
MBP field strengths at the photosphere and at the merging height are 1430 G and 120.4 G, respectively.^ 
The upper, monolithic part of the magnetic structure {z > Zm) is assumed to be cylindrically symmetric 
and is described by a modified version of the magnctostatic flux tube model of Hasan et al. (2003). In 
the present case, the flux tube is contained within a cylinder of transverse radius Sq = 10 Mm, which 
simulates the effect of neighboring network elements. The internal and external gas pressures are taken from 
semiempirical VAL/FAL models by Vernazza, Avrett, & Loeser (1981) and Fontenla, Avrett, & Loeser (1990, 
1991, 1993, 2002). The magnetic field components Bs{s, z) and Bz{s, z) in cylindrical geometry (using s as 
the perpendicular distance from the central flux tube axis) are computed by varying the shapes of the field 
lines until a minimum-energy state is obtained (for details see Hasan et al. 2003). The horizontal distribution 
of Bz{s,Zm) at the merging height is adjusted iteratively such that the magnitude of the field B{s,Zm) is 
independent of s. This is needed for consistency with the constant field strength of the thin flux tubes just 
below the merging height. 

Figure 2 shows the final iterated magnetic structure. The cross sections of the flux tubes below the 
merging height {z < 0.6 Mm) were computed by matching Bg and Bz at the merging height and assuming 
that Bg is independent of height within each flux tube (different tubes have different Bg). The transverse 
radius of the "network patch" in the photosphere is about 3 Mm. Wc can estimate the number N of MBPs 
inside the approximately circular patch, but we note that the wave analysis below does not depend on the 
value of N. From the conservation of magnetic flux, we know that the filling factor of MBP flux tubes 
at the photosphere is given by B{zm)/ Bq k. 0.08. Observationally, the transverse radii of MBPs at the 
photosphere are 50 to 100 km. The number of flux tubes required to fill an area that is 8% of the total patch 
area [A = ttSq , with s© = 3 Mm) ranges between N !v 75 (for the upper limit on MBP radius) and N « 300 
(for the lower limit). We can also define a mean distance between nearest neighbors rf„„ by dividing up up 
the network patch into N equal-area circles and defining dnn as twice their radii. Thus, rf„„ = 2sqN~^/'^ , 
which gives values of 350 and 700 km for the above limiting cases. 

Before continuing, it is worthwhile to ask the question: "Why is it not possible to extend the 'mono- 
lithic' fiux-tube model all the way down to the photosphere?" Observationally, the photosphere is far from 
monolithic in its distribution of magnetic fields. Granulation causes both the formation of thin MBP flux 
tubes and their spreading out (via random walk). The larger-scale supergranulation pushes the magnetic 
elements back together into network lanes and vertices. This competition between MBP spreading and con- 
vergence leads to a dynamical equilibrium described by the different assumptions applied above and below 
the merging height. A more complete model must contain physics that naturally captures this equilibrium 
state, but hopefully without the abrupt transition assumed here at Zm- 

3.2. Superradial Expansion in the Solar Wind 

Here we describe the plasma parameters along the central axis of the modeled network element at 
heights ranging from 12 Mm (the top of the magnetostatic model) to interplanetary space. We adopt a 



^ The magnetic field is assumed to point outward — i.e. , it is assumed to be of northern polarity — but the physics is the same 
for an inward pointing field. 
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slightly modified version of the polar magnetic field configuration of Banaszkiewicz et al. (1998): 

Bo _ 2 4.281 1 0.67 

1.789 G ~ x^^ aJxTocf ^ exp[384(a: - xq)] 



(2) 



where Bq is now defined as the radial component of the field along the axis of symmetry, x = t/Rq, and the 
above description applies only for x > xq, with xq = 1.0172 (i.e., the top height [12 Mm] of the magnetostatic 
grid). The model above uses the same current-sheet constant (oc = 1.538) used by Banaszkiewicz et al. 
(1998), but with a 5% modification to their preferred quadrupole constant. We also add an exponential 

correction term (which drops rapidly to nearly zero for x > 1.03) to ensure that the value and slope of Bq 
match those of the magnetostatic model at xq- Figure 3a plots the product x^Bq versus height both below 
and above a;o- This quantity is proportional to the inverse of the traditional superradial divergence factor /, 
and it shows the outer monopolar expansion region {Bq oc r^^) as a constant. For comparison we also plot 
the analytic functions used in the funnel models of Hackenberg, Marsch, & Mann (2000) and Li (2003). 

For the radial dependence of the electron density, we use a function motivated by fits to white-light 
polarization brightness measurements in the extended corona: 

rie _ 1 25 300 1500 5796 

~ ^ ^ ^8 ^16 ^33.9 V<^) 



1.3 X IQS cm-3 a;2 a;^ 

which applies only for x > xq- The overall scale set by the inverse square term was adjusted to match in situ 
density measurements at 1 AU (i.e., x = 215). The middle three terms above were adjusted to produce 

agreement with measurements by, e.g., Guhathakurta & Holzer (1994), Fisher & Guhathakurta (1995), 
Doyle, Teriaca, & Banerjee (1999), Esser & Sasselov (1999), and Figure 10 of Lie-Svendsen, Hansteen, & 
Leer (2003). The last term was set to match the value and slope of the magnetostatic model density at xq. 
Figure 3b plots the hydrogen number density nn, computed assuming a helium-to-hydrogen number density 
ratio of 0.05 (i.e., rig = 1-1 nn, with the total mass density p given by 1.2niimu)- 

Figure 3c shows several velocity quantities along the central axis of the fiux tube. The hydrogen outflow 
speed u was computed by mass fiux conservation (i.e., pu/Bo = constant), with the mass loss rate set by 
Ulysses measurements in interplanetary space (Goldstein et al. 1996). At 1 AU, the product unn = 2 x 10^ 
cm~^ s~^, and thus we compute u — 781.2 km s^^ at 1 AU. At an infinite distance, u approaches a constant 
value of 781.9 km s~^. The Alfven speed, defined as 

Va = BqI^p , (4) 

is nearly constant below the merging height (and would be precisely constant for an isothermal atmosphere) 
and rises to two successive maxima: 2530 km s""'^ at r = 1.004 ii©, and 2890 km at r = 1.53 i?0. The 
Alfven speed drops to 31.3 km s~^ at 1 AU and decreases nearly exactly as 1/r after that. The Alfven 
critical point, where u = Va, is at ta — 9.70 Rq. 



4. Photospheric Fluctuation Spectrum 

The lower boundary condition for our model of Alfvenic fiuctuations is the power spectrum of transverse 
MBP motions in the photosphere. The dynamical behavior of G-band bright points has been studied 
observationally by a number of groups (see, e.g., Muller et al. 1994; Berger & Title 1996; van Ballegooijen 
et al. 1998; Berger et al. 1998; Krishnakumar & Venkatakrishnan 1999; Nisenson et al. 2003). In this section 
we present an empirical description of MBP dynamics as a linear superposition of two types of motion: (1) 
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the "random walk" undertaken by isolated flux tubes, and (2) a series of rapid "jumps" that occur when 
individual flux tubes merge, fragment, or reconnect with surrounding magnetic field. 

The general procedure for specifying the power spectrum of horizontal MBP kinetic energy (as a function 
of frequency) is illustrated in Figure 4. The primary measured quantity is a time series of discrete position 
measurements for the MBPs which can be differentiated to obtain horizontal velocity components {t) and 
Vy{t), with the z direction being normal to the solar surface. (Observations support the rational assumption 
that there is no preferred global direction on granular scales, so below we discuss just the component and 
assume that Vy is statistically equivalent.) The timescale dependence of MBP fluctuations is encapsulated 
in the velocity autocorrelation function, which we define for a time sequence of A'' measured velocities w^.n 
as 

^ N-m 

^xx,m ~ T7 / . '^x.n '^x.n+m (5) 

JN — m ^ — ' 

n— 1 

for the time index n between 1 and N and an arbitrary delay time represented by index m. This expression 
is the discrete version of the more general definition 



Cxx{t) = lim — / dtv^{t)vx{t + T) (6) 



for delay time r. The unidirectional power spectrum (i.e., \vx\'^ per unit interval of frequency uj) is the 
Fourier transform of the autocorrelation function, with 

Px{u) = y^°° dT C,,(r) e'^^ (7) 

(e.g., the Wiener-Khinchin theorem; see also van Ballegooijcn et al. 1998). Note that Px is defined for both 

positive and negative w. For the present applications, all functions are symmetric about zero frequency, and 
we will thus consider only positive frequencies (simultancoiisly doubling the normalization of Px to conserve 
total energy). In general, we specify the kinetic energy power spectrum Pk, which is defined in such a way 
as to integrate to the total kinetic energy density Uk in transverse motions: 



Uk = — 



with 



/ duPKiuj) , (8) 
Jo 

Pk{uj) = I [2Pxiuj) + 2Py{uj)] = 2pPx{uj) . (9) 

The factors of two inside the square brackets take account of the negative frequencies, and the last expression 
assumes Px = Py 

In the analysis below we derive power spectrum components for the two assumed phases of the MBP 
motion: the random walk (subscript w) of isolated flux tubes, and the occasional discrete jumps (subscript 
j) caused by merging, fragmenting, or reconnecting. Specifying the power separately for these two phases is 
not the ideal solution, but it is all that can be done at present. Ideally, the proper observational procedure 
would bo to determine the com,plete time scries I'xif) for the walk and jump phases taken together, then 
compute the autocorrelation function and total power spectrum consistently. Unfortunately, because MBPs 
fragment or merge during the jump phases, it is extremely difficult to "follow" a single feature during these 
times to determine the complete time series. 
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The random-walk component of MBPs, during the times they exist as separate entities, was studied 
by van Ballegooijen et al. (1998) and Nisenson et al. (2003). These observations yielded the result that the 
discretely derived autocorrelation functions can be fit well by Lorentzian functions (see Figure 4e), with 

and the more precise observations of Nisenson et al. (2003) gave values of cr^ 0.8 km^ s~^ and Tyj ~ 60 
s. We adopt these values in the models below. Thus, the walk-component of the kinetic energy spectrum is 
given by 

PkA'^) = polr^e-^-'^ . (11) 

This is plotted in Figure 4g, but note that all power spectra arc plotted as the product cjP(cj) because this 
denotes the power per decade of frequency (i.e., the energy density per unit logw). Maxima in this quantity 
highlight the frequencies that contribute most to the total wave energy. 

The impulsive "jump" phase of MBP motions is described by Choudhuri, Dikpati, & Banerjee (1993), 
Berger et al. (1998), Hasan, Kalkofen, & van Ballegooijen (2000), and others. The potential for rapid 

transitions in the locations of thin flux tubes is also indicated in empirical models of the quasi-equilibrium 
evolution of granular magnetic fields (van Ballegooijen & Hasan 2003), in which slow motions of separatrix 
surfaces in the photosphere are amplified at larger heights due to flux-tube expansion. We model an impulsive 
MBP event as narrow Gaussian enhancement in Vxit) centered on an arbitrary t = 0. A series of these 
events is assumed to occur with a mean time interval Ai between events.'^ The Gaussian jump, with velocity 
amplitude CTj and 1/e half- life Tj , is described by 

v^) = a,-e-(*/-')' (12) 

between times and -|-Af/2, and the limit ^ oo in eq. (6) is not taken. The autocorrelation 

function is thus 

^ 



(see Figure 4f) and the resulting kinetic energy spectrum is given by 



^2^2 / , ,2^2 



At 



PkA^) = ^exp -— ^ (14) 



(sec Figure 41i). In the models below we adopt At = 360 s and tj = 20 s, which are consistent with the 
observations of Berger et al. (1998). The velocity amplitude of the jump aj is known with much less certainty 
because it represents the "tail" of the observed distribution of speeds. Berger & Title (1996) found speeds 
up to 5 km s~^, so this seems to be a rough upper limit to aj. This is our only true free parameter. 

Because the walk and jump phases of MBP motion seem to be statistically uncorrelated, we compute 
the full kinetic energy spectrum Pk as the simple sum of Pk,w and Pkj (see Figure 4i). It is then possible 



^Although we believe a constant interval At captures the essential nature of these jumps and their contribution to the energy 
spectrum, a more accurate way of modeling them would be to sample from an empirically derived "waiting-time" probability 
distribution. 
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to integrate this spectrum over frequency (see eq. [8]) to obtain the transverse velocity variance of MBP 
motions in the photosphere: 



ml 



2al + \/2^a|r,/At 



1.6 + 0.139 



1 kms~i 



km^-2 



(15) 
(16) 



where the latter expression uses the values adopted above for a^, tj, and At. Note that even for a large 
impulsive velocity of aj ~ 6 km s~^, the root-mean-squared velocity is significantly smaller than this value 
{{SV)q « 2.5 km s~^) because the jumps occur infrequently. 

Finally, we also compute the total energy spectrum Ptot (i-e., including contributions from both kinetic 
and magnetic energy) at the photosphere using the analytic relations for isothermal thin flux tubes given in 
Appendix A (specifically, cq. [A8]). For very high frequencies the kinetic and magnetic energy components are 
in equipartition, and Ptot is just twice Pk- For very low frequencies the kink-mode waves are evanescent and 
the physically realistic solution contains much more kinetic energy than magnetic energy (thus, Ptot ~ Pk)- 
The resulting total energy spectrum, plotted as a solid line in Figure 4i, is essentially our lower boundary 
condition on the amplitudes of Alfven waves of various frequencies. We describe how this information is 
folded into the global solutions in § 5.3. 



5. Non-WKB Wave Analysis 

We model the transverse wave properties in the open magnetic regions described in § 3 as purely linear 
perturbations to the assumed zero-order background plasma state. The basic MHD equations that need to 
be solved are the mass and momentum conservation equations and the magnetic induction equation, given 

by 

§+V-(pv)=0 (17) 

^ + (vV)v = -ivp + g+-^[(VxB)xB] (18) 
at p Airp 

— =Vx(vxB) (19) 

where the velocity v and magnetic field B are not yet separated into zero-order and first-order parts, p is 
the gas pressure, and g is the gravitational acceleration. 

We apply these equations in two regions with very different physics: 

1. Below the merging height (§ 5.1) wc model the waves as incompressible Lagrangian perturbations of 
the central axis of a thin, strong-field flux tube that expands superradially and is surrounded by a 
field-free region. In this region we assume the outflow speed of the solar wind is negligibly small. In 
addition to the properties within the thin tube, we also specify the density external to the tube Pe in 
the fleld-free region. 

2. Above the merging height (§ 5.2) wc model the waves as incompressible Eulerian perturbations filling 
the volume of the expanding flux tube, which is assumed to be surrounded by similar tubes. All 
background properties, including a nonzero solar wind speed, vary only in the radial direction. 
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For both regions we transform the MHD equations above into wave equations with an assumed e"^* time 
dependence. The equations are then solved "monochromaticany" for a grid of frequencies, also assuming 
that in each solution the frequency remains constant as a function of height. We do not make the WKB (i.e., 
eikonal) approximation^ that wavelengths are small compared to background scale lengths; indeed we do not 

even need to define the concept of wavelength because the complete spatial oscillation pattern is computed 
numerically as a function of height. The solutions for individual frequencies are subsequently assembled 
into a full radially varying power spectrum, normalized by the empirically derived power spectrum at the 
photosphere (§ 4). 

The following subsections present the specific equations that we solve in the two regions outlined above. 



5.1. Thin Flux Tubes 

Below the merging height, the MBP flux tubes are shaken transversely and kink- mode waves are excited 
(see also Wilson 1979; Spruit 1981; Ulmschneider et al. 1991). For incompressible perturbations about the 
equilibrium state, the density p is a zero-order quantity, the velocity v is horizontal and a first-order quan- 
tity, and the magnetic field B has a zero-order vertical component and a first-order horizontal component. 
Following the motion of the thin tubes, we write the Lagrangian forms of the momentum and induction 
equations as follows: 

rfv / (B-V)B , , 

f = (B.V)v (21) 

where the advective derivative 

follows the motion of the tube's central axis. Above we have assumed that V • v = as a statement of mass 
conservation for incompressible flows. 

We write the scalar horizontal perturbations in velocity and magnetic field as v± and B±, and we 
implicitly assume linear polarization of the waves in a single transverse dimension. In many ways, though, 
the equations to be derived are degenerate with toroidal Alfven waves at larger heights (e.g., Heinemann 
& Olbert 1980) and our assumption should not appreciably limit the generality of the results. (Other 
polarization modes have been studied by, e.g., Spruit 1982, 1984; Lou 1993; Roberts 2000; Noble, Musielak, 
& Ulmschneider 2003; Ruderman 2003.) 

In order to derive the wave equations we note three aspects of thin tubes in the solar atmosphere: 

1. For an oscillating flux tube, the direction perpendicular to its instantaneous axis will, generally, be 
inclined with respect to the radial direction away from the Sun. The components of the above equations 
parallel to the flux tube axis are uninteresting, and will be considered to be "solved" by the zero-order 
background state described in § 3. 



*By "WKB" wc refer broadly to the use of an asymptotic expansion that facilitates the solution of the linear differential 
equations. Specifically, for the wave equations presented in this paper, the "WKB limit" is that of pure outward propagation 
with no reflection. The use of this acronym that cites the contributions of Wentzel, Kramers, and Brillouin does not imply 
neglect of the earlier work of Liouville, Green, Carlini, Rayleigh, Jeffries, and others. 
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2. We assume transverse total pressure balance between the flux tube and the surrounding field-free 
region, and that the field-free region is in simple hydrostatic equilibrium. Thus, 

^(p+^^ = "^Pe = Peg . (23) 

3. The motion of the tube induces motions in the surrounding field-free region, which in turn must have 
a back-reaction on the tube's original motion. Spruit (1981) took this into account by increasing the 
apparent inertia of the tube. Thus, for the perpendicular component of eq. (20), the factor of p in the 
advection term on the left-hand side must be replaced by (p+pe)- This is equivalent to the assumption 
that the tube "carries along" a parcel of the surrounding fluid with equal kinetic energy density to that 
of the tube itself. Osin, Volin, & Ulmschneider (1999) reviewed different approaches to the inclusion 
of this back-reaction effect and found that Spruit's (1981) approach adequately describes the physics 
for transverse oscillations of a nearly vertical flux tube. 



With the above considerations. Spruit (1981) showed that the perpendicular component of eq. (20) can 
be expressed as a linearized wave equation 

aV ^ gAp dv± 2 /24^ 
df^ Ptot dr dr^ ^ ' 

where the magnitude of the gravitational acceleration g = GMq/t^ is nearly constant over the heights we 
consider. In Spruit's (1981) ideal limit that the flux tubes are completely isolated from one another, the 
density quantities introduced above are defined as 

Ptot = P + Pe (25) 

Ap = p-pe (26) 

and V^h = Bo/^/inpu^ is a modified kink- mode phase speed that takes the surrounding inertia into account. 
Note that in the Lagrangian picture, v± is the time derivative of the horizontal displacement ^. The first 
term on the right-hand-side of eq. (24) is due to the buoyancy of the low-density flux tube. The second term 
on the right-hand-side is the magnetic tension restoring force due to the curvature of the flux tube. The 
Lagrangian induction equation is given simply by 

-w-^'^- ^^^^ 



Analytic solutions to the above equations are possible when the radial derivative terms have constant 
coefficients; this occurs in an exponential isothermal atmosphere (see Appendix A for details). Traditionally, 
the solutions to these equations display evanescence for frequencies below a critical cutoff value ojc- For the 
adopted background state at the photosphere {z = 0), pe/p = 2.35 and Vph = 6.672 km s~^, and Appendix A 
gives an analytic estimate for the corresponding critical period (2'7r/wc) of about 12.5 minutes. This period is 
significantly longer than the acoustic cutoff frequency of 3-6 minutes in the photosphere and chromosphere. 
Thus, the kink mode has been suspected for several decades as being able to transport more convective wave 
energy up to the corona than acoustic waves. 

Before discussing our numerical solutions to eqs. (24)-(27) between z = and z^,, one simplification 
assumed above must be reexamined. Just below the merging height, the flux tubes cannot be considered truly 
isolated from one another. The enhanced inertia assumed by Spruit (1981) assumed that the surrounding 
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fluid carried along by a given tube is all field-free, but near the merging height this is not the case. Spruit 
(1982) gave the equations for selected kink-mode wave properties for the general case where the surrounding 
medium has a nonzero field strength, but here we deal with the encroachment of neighboring flux tubes in 
a simpler manner. In the above equations we express ptot and Ap by 

Ptot = P+(l-0)Pe (28) 

Ap = (p-pe)(l-0) (29) 

whore (f> is essentially a statistical filling factor of neighbor tubes within the near-tube region that gets carried 
along with a tube's oscillatory motion. The isolated tube limit is = 0, but at the merging height (and 
above), (j) =1. The specific form of the modifications above were constrained by the need for both V^h = Va 
and for the buoyancy term in eq. (24) to vanish in the "merged" limit of </> = 1. The reduction of Ap also 
reduces the critical frequency for evanescence (see eq. [A2]), and at the merging height ujc 0. We derive 
0(z) by using magnetic flux conservation together with the assumption that the overall area subtended by 
the full network patch is constant between the photosphere and the merging height. Thus, 

^ { B,{z^)/B,{z) , z<Zm .30. 

11, Z>Zm 

where B(j{zm) = 120.4 G. Because Eq increases rapidly with decreasing height below Zm, <!> rapidly drops 
from 1 (at Zm = 600 km) to 0.5 at ^; = 450 km, then more slowly down to 0.08 at ^; = 0. In Figure 3 we plot 
the resulting values of ptot (as a hydrogen number density) and Vph below the merging height. 

The wave equation (eq. [24]) is simplified by assuming an e'"* time dependence with a real frequency 
u). We solve numerically for the radial dependence of v± by expressing the second-order wave equation as 
two coupled first-order ordinary differential equations in v± and dv±/dr (with both quantities assumed to 
be complex), and using fourth-order Runge-Kutta integration (e.g., Press et al. 1992). The upper boundary 
conditions at z^ are specified by the solutions of the wave equations above z„i (sec next section), and the 
numerical integration proceeds from Zm down to the photosphere. The coupling of solutions below and above 
the merging height is discussed in § 5.3. 



5.2. Wave Reflection in the Solar Wind 

Above the merging height, the transverse incompressible fiuctiiations act as MHD Alfven waves, and our 
solution procedure largely follows that of Heinemann & Olbert (1980) and Barkhudarov (1991). Formally, 
the mechanisms of WKB theory can be extended to describe linear wave reflection (e.g., HoUweg 1990), but 
we follow the usual non-WKB formalism in order to solve for the radial dependence of the transmitted and 
reflected wave properties from the merging height all the way to a distance of 4 AU. 

The monochromatic non-WKB wave transport equations are derived from the mass, momentum, and 
induction equations listed above (eqs. [17]-[19]) in the limits that all background quantities vary only in radius 
and that the velocity and magnetic field perturbations are perpendicular to the zero-order field direction. 
Khabibrakhmanov & Summers (1997) showed how to treat general vector operations in a superradially 
expanding flux tube. We express the wave properties in terms of Elsasser (1950) variables, deflned here as 



(31) 
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(see also Tu & Marsch 1995), with z_ representing outward propagating waves and z+ representing inward 
propagating waves. In terms of these variables, the incompressible first-order equations are expressed as two 
coupled transport equations: 

where the (signed) scale heights are Hd = p/{dp/dr) and Ha = VA/{dVA/dr). These equations are valid 
for superradial divergence and for all values of the zero-order outflow speed u. Note that linear reflection 
arises because of the Zzf. term on the right-hand side. 

In Appendix B wc discuss additional details about these equations and how they are equivalent to other 
versions given in earlier work. To our knowledge, the "compact" form of cq. (32) has not been recognized 
fully, although the expressions of Heinemann & Olbert (1980) and Khabibrakhmanov & Summers (1997) 
were closely related. The above form of eq. (32) is particularly useful in showing the large-scale density 
dependence of the wave amplitude in two limiting cases when the outward propagating waves are dominant 
(i.e., \z-\ » Near the Sun, where u -C Va, we can approximate the lower-sign version of the above 

equation as 

dr Ap dr ^ ^ ^ 

and thus is proportional to p^^/^ as predicted by WKB theory (see also Moran 2001). Similarly, far from 
the Sun, where u ^ Va, z_ is seen to be proportional to More general properties of the non-WKB 

solutions of eq. (32) are discussed in § 6, and also by MacGregor & Charbonneau (1994) and Krogulec et al. 
(1994). 

As with the solutions below the merging height, we assume an oscillatory time dependence of the Elsasser 
variables (i.e., e*"* with a real, constant frequency), and wc solve for their radial dependence numerically. 
We restrict our solutions to positive frequencies and note that taking the negative of a given frequency 
produces solutions to eq. (32) that are the complex conjugates of the analogous solutions obtained with 
w > 0. Thus, the radial evolution of physical quantities (i.e., real wave amplitudes) is unaffected by the 
sign of w. The existence of the Alfven critical point complicates the numerical solution of eq. (32), but we 
follow the general solution procedure outlined by Barkhudarov (1991). Once the oscillatory time dependence 
has been assumed, the complex amplitudes z+ and Z- are expressed as the products of real amplitudes and 
phase factors of unit magnitude. There arc then four ordinary differential equations for these quantities 
that are solved first at the Alfven critical point ta (analytically, using certain physicality constraints such 
as the requirement that the outward wave energy always exceed the inward wave energy), then we integrate 
numerically using the fourth-order Runge-Kutta method both upwards and downwards from ta- Some of 
Barkhudarov's (1991) expressions had to be modified to take account of the superradial divergence of the 
magnetic field. The linear amplitudes \z±\ are specified only to within an arbitrary normalization factor, 
although both the phase factors and all relative quantities (such as ratios of Elsasser amplitudes at different 
radii) do not depend on this normalization. 



5.3. Solution of the Coupled Wave Equations 

Our baseline model consists of a grid of 300 frequencies, evenly spaced in logo; over five orders of 
magnitude with periods ranging from 3 to 300,000 s. The discrete radial grid extends from just above 
the photosphere [z = 10 km) to 4 AU {z = 859 Rq) and contains 11809 grid points distributed mainly 
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logarithmically, but with some regions (like the transition region) sampled more finely. In the photosphere, 
chromosphere, and low corona (below xo), the relative grid separation Az/z is 0.00064. Within the most 
rapidly changing 20 km of the transition region, Az/z is decreased to 0.0001. In the extended corona and 
solar wind, Az/z is made to gradually increase to 0.016 at the outer boundary. The Runge-Kutta algorithm 

also has an adaptive stcpsizc that subdivides the above grid zones until a relative accuracy of 10~^ is achieved 
in the integration variables. This degree of accuracy is needed to follow the oscillatory behavior of the waves. 

The non-WKB wave equations are solved first for each frequency in the upper corona/wind region as 
described in § 5.2, and the resulting Elsasser variables at the lower boundary (i.e., the merging height) are 
used to compute the complex values of v± and B± at that height. These quantities, still with an arbitrary 
degree of normalization, are used as the upper boundary conditions for the numerical solution of the flux- 
tube wave equations given in § 5.1. The induction equation, eq. (27), is used only to convert the boundary 
condition for B± into a condition for dv±/dr. We make the assumption that all of the Alfvenic wave 
energy in the upper region is converted smoothly into kink-mode wave energy in the lower region. After the 
transport equations arc solved in both regions, the photosphcric MBP power spectrum (derived in § 4) is 
used to renormalize the wave power quantities at all heights. To show how this is done we must first define 
the kinetic and magnetic energy densities for each monochromatic model: 

Ek = , Eb = (34) 

and also the supplementary quantities 

E^^'-^, E^^e^^^^=E^-E_ (35) 



where below the merging height we use ptot for the above densities. The energy densities defined above do 
not depend on time t since they contain products of e**^* and its complex conjugate. These definitions also 
ensure that the total fluctuation energy for each frequency satisfies £^tot = Ek + Eb = E+ + E-. In the 
simple WKB theory (i.e., all outward propagating waves), E+ = 0, Ek = Eb, and Etot = E- = —Ex, and 
the departure from WKB theory can be assessed roughly by departures from these ideal energy identities. 

In § 4 we derived the total power Ptot(w, Rq) at the photosphere. The numerical integrations described 
above gave us the various E quantities on a two-dimensional grid in oj and radius (with the energies for 
each frequency known up to an arbitrary multiplicative constant). We thus compute "renormalized" power 
spectra on the discrete grid as 



PsiuJ.r) I ^ Ptot{Lu,RQ) ^ I EBiLU,r) 
P±ito,r) I Etot{oj,RQ) I E±{LU,r) 
ftot(w,r) J [ Etoti(^,r) , 



(36) 



Once the power spectra have been defined we can then specify various frequency-averaged energy densities. 
In general, we define 

/•CxD 

C/tot = / dwPtot(c^) (37) 



^0 

with analogous definitions for Uk, Ub, and U± (see also eq. [8]). For ease of interpretation we also define the 
frequency-averaged velocity, magnetic, and Elsasser amplitudes as the quantities in angle brackets below: 

^, = £M!, = (38) 



(39) 



6. Results and Observational Implications 

6.1. Linear Wave Properties 

The procedures outlined above resulted in a large amount of numerical data (~350 megabytes) describing 
the behavior of non-WKB kink-mode and Alfven waves as a function of frequency and radius. In this 
section we attempt to distill and present the salient results in three gradiial steps: (1) in Figures 5 7 we 
present frequency-dependent wave properties that have not yet been renormalized to the photospheric power 
spectrum (and thus are plotted as dimensionless ratios), (2) in Figure 8 we show the total power spectrum 
as a function of frequency for selected radii, and (3) in Figures 9-18 we present various frequency-integrated 
quantities that depend on the photospheric power normalization. 

In order to determine to what degree the waves in various regions depart from ideal WKB theory, we 
show in Figure 5 the outward propagating wave action flux p (defined in Appendix B) for a selection of 
periods, each normalized to the value of p in the photosphere. For fluctuations having wavelengths so 
small that the local plasma parameters arc approximately constant over several wavelengths, wave action is 
conserved and the quantity / is constant. Note from Figure 5 that in the extended corona and solar wind 
(e.g., ^; > 1 Rq) Alfven waves with periods shorter than a few hours obey wave action conservation, but for 
periods exceeding -^10 hours this breaks down. Note that for all computed periods the wave action is not 
conserved below the transition region [z w 0.003 i?©); this region acts as a sufficiently sharp "barrier" to 
induce significant reflection in the chromosphere (see also Wcntzel 1978; HoUwcg 1981; Campos & Gil 1999). 

In much of the previous work on non-WKB Alfven wave reflection, the departures from WKB theory 
have been characterized as a function of frequency. For frequencies exceeding a critical value uja-, the 
resulting wavelengths are so short that the wave propagates as if it were in a homogeneous medium and 

there is negligible reflection. For frequencies lower than loa there is significant refiection, and as a; ^ 
the oscillation approaches the properties of a standing wave with equal inward and outward power. In a 
magnetized hydrostatic atmosphere, the critical frequency is given by the local value of |9VA/9r|, or more 
precisely for arbitrary expansion factors, |V • V^l (e.g., Ferraro & Plumpton 1958; An et al. 1990). In a 
supersonic wind, though, the radial dependence of the Alfven speed is no longer the dominant factor in 
determining how much reflection takes place. At large distances from the Sun, Heinemann & Olbert (1980) 
and Barkhudarov (1991) showed that u)a is given approximately by Uoo/2rA, where Uoo is the asymptotic 
outflow speed and ta is the Alfven radius. (We use this expression as a fiducial definition of w^i, but we also 
note that it neglects several order- unity correction factors that depend on the fiow tube geometry.) For the 
zero-order solar wind model defined in § 3.2 we find an associated critical period 2-^ /uja of about 30 hours, 
or 1.25 days. 

Figure 6 shows contours of the Alfven ratio — i.e., the ratio of kinetic to magnetic energy density Ek/Eb 
as a function of period and height. For ideal MHD Alfven waves in a homogeneous medium this ratio is 1 
and the waves are in energy cquipartition (Walen 1944). Figure 6 can be broken up into four "quadrants" 
that have the following limiting properties: 

1. Upper left: For long periods below the transition region, most of the wave energy is kinetic with only 
a negligible magnetic energy density. This is consistent with the predictions of kink-mode wave theory 
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for the so-called "shallow" evanescent solution. In Appendix A we discuss several reasons why the 
solar atmosphere is suspected to naturally prefer this solution. 

2. Lower left,: For short periods below the transition region, the Alfvcn ratio rapidly fluctuates above and 
below 1. This occurs because the wavelengths are small compared to the photospheric and chromo- 
spheric scale heights, but they are large compared to the scale heights in the transition region. Thus 
significant reflection occurred and there is a superposition of upward and downward waves. The fluc- 
tuations are well described by standing waves having flxcd nodes in v±^ and B±^ (HoUwcg 1981, 1984). 
If one could separate the solutions into the component upward and downward waves, their individual 
Alfven ratios would be -^1 as predicted by the kink- mode theory in Appendix A. (This can also be 
seen roughly by averaging over several nodes.) 

3. Upper right: For long periods above the transition region, most of the wave energy is magnetic as was 

found by other non-WKB solar wind models (e.g., Hcincmann & Olbert 1980). These long periods 
correspond to "quasi-static" motions of the field lines and have a similarity to the motions invoked in 
DC theories of coronal heating (e.g., Kuperus et al. 1981; van Ballcgooijen 1986; Milano, Gomez, & 
Martens 1997). 

4. Lower right: For short periods above the transition region, the plasma appears homogeneous to the 
relatively small- wavelength fluctuations and the ideal MHD equipartition holds. 



Figure 7 shows contours of the ratio Ej^/E_ which can be thought of as an effective local reflection 
coefficient. The ideal WKB theory corresponds to a ratio of zero (i.e., all outward propagation). There 
is significant reflection below the transition region, with a maximum value of this ratio of 0.99825 at the 
photosphere for a period of ~40 minutes. As discussed in § 5.2, we chose regularity conditions at the Alfven 
critical point that ensured more outward than inward power at all heights, and thus E+/E- must always 
be less than one. Below the transition region, the behavior of E+/E- versus frequency resembles analytic 
solutions that take account of the chromospheric stratification and the strong reflection at the transition 
region (e.g., HoUweg 1978a; Schwartz, Cally, & Bel 1984). Because of the finite size of the atmosphere below 
the transition region, a mild resonance structure is seen in the amount of reflected wave flux (not apparent 
in Figure 7). The resonances are not as sharp as in the isothermal models of Schwartz et al. (1984), though, 
because the nonzero temperature gradient results in a "smearing" of the preferred resonance frequencies. 

Above the transition region, the solutions plotted in Figure 7 arc very similar to those of Heinemann & 
Olbert (1980) and Barkhudarov (1991); there is significant refiection when cu < uja and much less reflection 
for higher frequencies (lower periods). For frequencies above to a — which we see below encompasses most of 
the dominant part of the power spectrum — we can use the analytic regularity conditions presented in § 4 of 
Barkhudarov (1991) to estimate the ratio E+/E- at the Alfven critical point: 



E+{uj,rA) _ J_ 
E-{u;,rA) ~ 



dVA 



dr 



for w > • (40) 



By comparing to the numerical results in Figure 7 we verify that this relation is accurate for periods less 
than ^200 minutes. It may be possible to also estimate the radial dependence of this ratio using similar 
analytic formulae. We defer this to future work, but note that this would be useful to models of coronal 
heating via turbulent cascade (see § 6.2.2). 

Figure 8 shows the result of convolving the above results with the photospheric power spectrum derived 
in § 4. In this figure we plot the total power Ptot('^) for selected heights, and with a specific choice for the 
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parainetc'i' (Tj of 3 km s"'^ (see below). The amount of power loss between the photosphere and the merging 
height is consistent with the analytic predictions of isothermal kink-mode theory. For high frequencies 
(i.e., periods less than 3-5 minutes), Ptot decreases from 2; = to by a factor of 15 to 20; this is 
approximately the level of the decrease of the background magnetic field between those heights, implying 
that the propagating kink-mode relation pv\ oc Bq applies. For low frequencies (presumably evanescent), 
Ptot decreases by a factor of 190, which is very nearly equal to the drop in density from z = Q to Zm- 
This is consistent with the prediction of nearly constant v±_ for the shallow evanescent solution discussed in 
Appendix A (and thus pv^ oc p). The somewhat irregular structure that develops in the spectrum between 
the photosphere and the merging height (and is passively advected outwards above z„i) is not numerical 
noise. Because of the nonisothermal temperature structure of the low chromosphere and the use of a radially 
varying filling factor (j), the properties of the waves below the merging height depend on frequency in a more 
complicated way than in the ideal case described in Appendix A. If the upward and downward propagating 
waves had identical strengths and dispersive properties, the standing- wave nodes in Ej^ and Eb would cancel 
out exactly in £^tot- However, ii^tot exhibits a weak nodal structure because of incomplete cancellation and 
thus contributes to the irregularity in the power spectrum (see also Schwartz et al. 1984). 

As in Figure 4, we plot in Figure 8 the product wPtot versus wave period to more clearly show the 
periods that make the greatest contribution to the overall wave energy. We quantify this concept by defining 
the averaged, or first-moment frequency as 



We plot this quantity versus height in Figure 14 below, but here we just note that the effective period 27r/ (w) 
is about 3.5 minutes in the photosphere, then it decreases to about 1.8 minutes above the merging height. 
Periods of a few minutes are natural to expect in the photosphere and chromosphere, and possibly in the 
corona as well (e.g., Chashei et al. 1999). However, it is reasonable to ask if these periods are expected to 
dominate in interplanetary space. In situ spacecraft generally measure fiuctuations in velocity, density, and 
the magnetic field with the most power at periods of a few hours (e.g., Goldstein et al. 1995). However, 
a spacecraft sitting still in the ecliptic plane at 1 AU would sec our model network flux tube rotate past 
in about 4 to 5 hours, with its component flux tubes (each originating in a different MBP) rotating by in 
substantially less than one hour. Thus it is possible that if the in situ power spectra actually do sample a 
"fossil" spectrum from the Sun, the dominant periods of order 1 hour could be due to the passage of many 
uncorrclatcd flux bundles past the spacecraft, and not the waves within any one flux bundle. (The ideal test 
would be to see if a spacecraft corotating with the solar rotation period measures a significantly different 
fluctuation spectrum than has already been observed.) 

Figure 9, a key condensation of results for this paper, plots the frequency-integrated velocity amplitude 
{5V) as a function of height for a selection of aj parameters and compares it to several different measurements 
of wave amplitudes. We discuss each set of measurements briefly, by number, below: 

1. The dotted line shows a best-fit height dependence for the microturbulence needed to match pho- 
tospheric and chromospheric line widths in the semiempirical VAL/FAL models (E. Avrett, 2003, 
personal communication; see also Fontenla et al. 1993, 2002). 

2. The filled circles show similar "nonthermal" line-broadening velocities measured (on the solar disk) 
in the transition region and low corona by the SUMER (Solar Ultraviolet Measurements of Emitted 
Radiation; Wilhelm et al. 1995) instrument on SOHO (the Solar and Heliospheric Observatory). The 
height of each point has been estimated by matching the height-dependence of temperature in the 
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VAL/FAL model with the assumed formation temperatures of the ions that correspond to each point 
(Chae et al. 1998). 

3. The crosses show nonthermal velocities inferred by SUMER measurements made above the solar limb 
(Banerjee et al. 1998). Off-limb Doppler broadening observations are better suited for measuring the 
properties of transverse Alfven waves than observations on the solar disk. 

4. The gray region shows lower and upper limits on the nonthermal velocity (Esser et al. 1999) as computed 
from off-limb measurements made by the UVCS (Ultraviolet Coronagraph Spectrometer) instrument 
on SOHO (Kohl et al. 1995, 1997). 

5. The stars show early measurements (Armstrong & Woo 1981) of the random wavelike component of 

the solar wind velocity from interplanetary scintillation observations of radio signals passing through 
the corona, with the advecting diffraction pattern being measured by more than one receiver. 

6. The error bars show a more recent determination of velocity fluctuations — specifically transverse to 
the radial direction — from radio scintillations (Canals et al. 2002) in the fast solar wind. 

7. The Helios and Ulysses probes measured time-averaged Elsasser amplitudes that we converted to a 
representative velocity amplitude using eqs. (38)-(39) and assuming Uk = i/s in the heliosphere. The 
data from both spacecraft were summarized by Bavassano, Pietropaolo, & Bruno (2000). 

Measured line-of-sight (one-dimensional) velocities have been multiplied by to take account of equivalent 
fluctuations in both transverse dimensions (e.g., eq. [9]). 

Because measurements (1) and (2) above refer to motions mainly in the radial direction, we do not 
expect them to correspond to transverse kink-mode or Alfven waves; they are plotted mainly for heuristic 
comparison. Note, though, that in Figure 9 we also plot (for the (jj = 3 km s~^ case) a dashed curve that 
shows the frequency-averaged amplitude of the magnetic fluctuations in velocity units: 



which for ideal MHD equipartition would be exactly equal to {5V). It is possibly a coincidence that this 
quantity so closely matches the Chae et al. (1998) observations, but this agreement may contain information 
about the mode coupling between transverse and longitudinal waves in the transition region. 

The three choices for the MBP-jump velocity amplitude dj used in Figure 9 are 0, 3, and 6 km s~^. The 
middle value seems to best match the ofi-limb nonthermal line broadening measurements, and we will use 
this as a baseline value in most subsequent calculations and plots. Note that the in situ measurements fall 
well below all of the reasonable choices for Oj (even the lower limit assuming no jumps whatsoever) and they 
exhibit a steeper gradient than the undamped models. The heliospheric "deficit" of wave power, compared 
to most prior assumptions about the wave power in the solar atmosphere, is well known (see, e.g., Roberts 
1989; Lou 1993; Mancuso & Spangler 1999). In § 6.2.2 wc propose a potential solution of these discrepancies 
based on a particular theory of turbulent wave damping in the solar wind. 

Figure 10 shows additional details about the height dependence of the frequency-integrated wave proper- 
ties. The ratio Uk/Ub shows the region of strong departure from energy equipartition around the transition 
region, and that kinetic energy exceeds magnetic energy by a slight amount even down to the photosphere. 
The dimensionless magnetic amplitude {5B) / Bq is less then 1 over most of the computed range of heights 
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(thus justifying the hnear approximation), but exceeds 1 above r w 19i?0. Nonhnear calculations (e.g., Lau 
& Siregar 1996) find that {5B)/Bo does not grow so high in the heliosphere and may saturate at values close 
to 1. 

Figure 11 shows the height dependence of the two Elsasser amplitudes. Above the transition region, the 
frequency-integrated outward propagating amplitude {Z-) is much larger than the inward amplitude (Z^). 
Below the transition region, the amplitudes approach a constant ratio {Z_) / (Z^) « 1.027. We also plot the 
height dependence of Z_ that would be expected for ideal WKB wave action conservation (i.e., / = constant) 
and we scale it to the computed value of (Z_) at 4 AU. The WKB value departs slightly from the non-WKB 
value between the transition region and a height of ~0.01 Rq and is substantially smaller below the transition 
region where there is strong inward wave power. The undamped linear curves all disagree markedly with 
the in situ measurements of the Elsasser amplitudes (see, however, § 6.2.2). 

Lastly, in Figure 12 we plot the frequency- integrated energy flux density F for the baseline (7j = 3 km 
s~^ case. Below the merging height this quantity applies to the energy flux carried upward by the MBP kink 
modes, and above the merging height it applies to the volume-filling Alfven waves. The wave flux is defined 
generally as 



which is based on the monochromatic definition given by Heinemann & Olbert (1980). For outward prop- 
agating waves obeying wave action conservation, one can write F = {7tot(l-5u + Va)- In the lowest layers 
of the atmosphere one can also ignore the factors above that depend on u and thus the net flux depends 
mainly on the difference between U- and f/+. Because of the strong reflection below the transition region 
there is strong cancellation; i.e., the magnitude of the purely outward wave flux VaU- exceeds the net flux 
by a factor of about 20. This means that ~95% of the kink-mode wave energy below the transition region 
is "trapped," with only 5% escaping. 

Analysis of a series of models having a range of aj values yields a potentially useful explicit expression 

for the photospheric wave flux in the MBPs. Recalling that the random-walk and impulsive parts of the 
power spectrum are assumed to be linearly independent (e.g., eq. [16]), we have flt the numerically computed 
photospheric fluxes with the following function: 



with less than 0.1% uncertainty in the flt (mainly due to roundoff error in the above expression having only 

three significant figures). Thus, for the reasonable range of aj between and 6 km s~^, Fq ranges between 
about 10^ and 10^ erg cm~^ s~^ (see also Musielak & Ulmschneider 2001). Note that for any given value 
of aj the jump- term in the above expression has a greater relative contribution to the total flux than the 
corresponding jump-term in eq. (16) has to the photospheric velocity variance. This occurs because the bulk 
of the power in the jump motions is at higher frequencies than the walk motions and is less siisccptiblc to 
evanescent decay. The flux is a frequency-integrated quantity that incorporates this aspect of the solutions, 
whereas the velocity variance does not. 

We also plot in Figure 12 several other flux quantities. The "granule-averaged" wave flux is given by 
the quantity (pF (for the flux tube filling factor </>, see eq. [30]), and it differs from F itself only below the 
merging height. The quantity (j)F spreads the wave energy out evenly within the modeled network patch, 
and thus is a more relevant quantity to compare with observations and predictions of wave fluxes that arise 
in the network from the convection. Note that the photospheric value of (j)F = 2.7 x 10^ erg cm~^ s~^ is of 
similar magnitude as the acoustic wave flux of 5 x 10^ erg cm~^ s~^ predicted by Musielak et al. (1994). We 



F = u{Uk + 2Ub) + Va{U_ - U+) , 
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also plot a "sTipcrgranule-averaged" flux 4>cF, where (f)c is a filling factor for the large-scale funnel/canopy 
magnetic structure shown in Figure 2. (We define 0c only below the top height [12 Mm] of the magnetostatic 
grid, as the ratio of the field strength at that top height [11.8 G] to the field strength along the central axis 
of the flux tube at lower heights.) The quantity (l)cF is what one would use in order to compute the total 
wave power (in erg s~^) integrated over the entire coronal hole. 



6.2. Wave Dissipation 

In this section we discuss two separate mechanisms that have been proposed to damp Alfven waves in the 
solar atmosphere and solar wind: viscosity (§ 6.2.1) and MHD turbulent cascade (§ 6.2.2). For clarity we do 
not treat these mechanisms together, though in a completely self-consistent model all damping mechanisms 
should be included and allowed to interact with one another. 



6.2.1. Linear Viscous Dissipation 

The flrst damping mechanisms proposed for MHD waves in the solar corona were collisional in nature 
(Alfven 1947; Osterbrock 1961). In the high-density plasma near the Sun, MHD waves can be damped 
by viscosity, thermal conductivity, electrical resistivity (i.e., Joule/Ohmic heating), and ion-neutral friction. 
For Alfven and fast-mode waves propagating parallel to the magnetic field, though, the dominant collisional 
dissipation channel in the corona is believed to be proton viscosity. Here we estimate the viscous damping 
expected for the background plasma conditions described above and find it to be negligible. 

We compute a collisional damping length Lc as a function of height z, and there should be appreciable 
damping only where Lc z. Generally, the product of the damping length and a linear damping rate 7 
(i.e., the imaginary part of the frequency) is the wave group velocity Vgr in the inertial frame, given here 
for dispersionless Alfven waves as w + Va. For viscous damping, 7 = z/pfc^, where Up is the proton kinematic 
viscosity and k is the radial wavenumber (van de Hulst 1951; Osterbrock 1961). Thus, for Alfven waves 
propagating upward and parallel to the background magnetic field, this leads to the general expression 
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(Tu 1984; Whang 1980, 1997). In a strongly collisional plasma the kinematic viscosity is given by the classical 
Braginskii (1965) formahsm (see also HoUweg 1986a). In the coUisionless solar wind, though, there is no clear 
prescription for the effective viscosity. We thus define the proton kinematic viscosity phenomenologically as 

Up = WpTeS = (46) 

TTlp 

where Wp is the proton most-probable speed, ks is Boltzmann's constant, Tp is the proton temperature, nip 
is the proton mass, and Tes is an effective viscous timescale. We adopt a fiducial proton temperature law 
for this analysis that comes from the semiempirical VAL/FAL model (used for the magnetostatic flux tube 
model) below x = 1.0172, and we use 

^ " 0.2 -h 0.02x0-8 0.21a;-33 ^ ' 

above 1.0172 Rq. This expression reasonably reproduces the proton temperatures measured by Helios and 
Ulysses in the highest-speed solar wind (Tu, Freeman & Lopez 1989; Goldstein et al. 1996) and it has a peak 
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in the extended corona of ^2.2 MK in general agreement with UVCS/ SOHO H I Lya measurements (e.g., 
Kohl et al. 1998; Esser et al. 1999). 

For strong coUisional coupling, the effective viscous time scale Tefj is given by the mean time between 
collisions (Braginskii 1965), 



where e is the proton/electron charge and In Ac is the Coulomb logarithm, taken here to be 21. The 
above timescale applies to the viscous damping of motions along the magnetic field. For shear motions 
transverse to the field, though, the appropriate viscous time (even in the limit of strong collisions) is reduced. 
The off-diagonal Braginskii coefficients in the stress tensor are obtained approximately by dividing TcoU by 
dimensionless factors that take account of the magnetic field. Thus, 



where m is either 1 or 2, and Hp is the proton Larmor frequency eBo/nipC. The w = 0, 1, 2 terms are roughly 
analogous to the direct, Hall, and Pedersen components of the MHD conductivity in an ionized plasma. 

In the low-density collisionless limit, the classical m = case cannot apply because it would imply the 
viscosity becomes infinitely large as the mean time between collisions becomes infinite (the "molasses" limit). 
Williams (1995) argued essentially that the general viscosity in a coUisional or collisionless plasma is found 
by taking the shorter of cither the m = or the m = 1 timescales. We assert, though, that the m = 2 case 
seems to be the most appropriate for the viscous damping of transverse Alfven waves in the collisionless 
solar wind. Consider the viscosity as an effective diffusion coefficient (i^/t) describing scattering events 
with mean-frce-path length and time scales of £ and t, respectively. The actual energy losses arise from the 
interparticle collisions with timescale t « Tcoh, but spatially the transverse structure would be dominated by 
features with sizes of order the proton thermal Larmor radius, and thus i ~ Wp/^lp. For these values the 
viscosity reproduces the m = 2 case. 

Figure 13 shows the radial dependence of the linear damping length for the three cases, m = 0, 1, 

and 2, computed for a wave period of 1 minute. For heights below about 0.3 Rq above the photosphere, 
it is unlikely that any damping would take place because all three cases have damping lengths much longer 
than the local height. This mid-corona distance (r « 1.3 to 2 Rq) is approximately where collisions start to 
become unimportant in coupling together electron, proton, and heavy ion plasma properties (e.g., Cranmer, 
Field, & Kohl 1999). Above this height, then, the appropriate damping length should transition to either 
the m = 1 or the m = 2 case, both of which are substantially longer than the local height in the corona. 
In the far solar wind, we believe the m = 2 case is the most realistic, and it always remains several orders 
of magnitude larger than the local height. Thus, our preliminary conclusion is that for both waves near 
the peak of the power spectrum (periods of a one to a few minutes), and for longer periods, linear viscous 
damping is unimportant as a significant attenuation mechanism for Alfven waves in the corona and fast solar 
wind. For periods much shorter than 1 minute, though, this kind of damping could be important — and may 
have already been responsible for the sharp drop in power inferred from photospheric MBP motions between 
0.1 and 1 minute. 
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6.2.2. Nonlinear Turbulent Dissipation 

The second type of wave damping we consider is turbulent dissipation; i.e., a nonlinear cascade of energy 
from large to small scales that terminates in an irreversible conversion of wave energy into heat. The need 
to include some kind of nonlinear damping or saturation can be seen in Figure 10, where above a distance of 
~20 Rq the magnetic fluctuation amplitude begins to exceed the background field strength (in opposition 
to in situ observations) and the linear assumption breaks down. Addressing the full problem of anisotropic 
MHD turbulence in the solar corona and solar wind, though, is well beyond the scope of this paper (see, 
though, Cranmer & van Ballegooijen 2003 for a summary of many related issues). Here we include only one 
specific aspect of this theory: a phenomenological damping rate that depends on the properties of the largest 
scales in the turbulence and not on the details of the cascade. 

The presence of both outward and inward propagating Alfven waves is a necessary prerequisite for the 
nonlinear coupling to higher wavenumber that drives the cascade (see, e.g., Matthaeus et al. 1999; Dmitruk, 
Milano, & Matthaeus 2001). Any expression for the turbulent energy decay must then depend on both the Z- 
and z+ Elsasser amplitudes, and thus the inclusion of turbulence introduces further coupling between these 
two modes. A straightforward phenomenological form for the nonlinear transport has been suggested by Zhou 
& Matthaeus (1990) from the standpoint of "reduced MHD" (RMHD; see also Strauss 1976; Montgomery 
1982). Effectively, the transport equation is the same as eq. (32), but with the following term added to the 
right-hand side: 

(50) 

where L± is a transverse length scale representing an effective correlation length of the turbulence — i.e., a 
similarity length scale associated with the energy transport from large to small eddies. (It may be imprecise 
to refer to this length scale as an "outer scale," but for simplicity we treat this as a synonym.) We assume 

— 1/2 

that L_L scales with the transverse width of the open flux tube; i.e., that it remains proportional to Bq 
Speciflcally, we normalize L± by specifying its value at the merging height Zm- We also define the constant 
A as the ratio of L± at the merging height to the transverse radius of the network flux bundle at that height: 
3 Mm (see Figure 2). Equivalently, we can express 

33 A 

L^{r) « -y== Mm (51) 

if Bq is measured in Gauss (see also Hollweg 1986b). 

In order to know how best to implement the nonlinear term given above, we need to understand over 
what range of heights this term acts to appreciably damp the waves. We do not expect a strong RMHD 
cascade to develop in the region below the merging height because — despite the stochastic nature of the 
fluctuation spectrum — the flux tubes do not strongly interact with one another. For heights above Zm we 
can compare the effective damping length L± to the local height z above the photosphere. For A = 1, we find 
that L± < z only above a height of z « 0.01 Rq. It makes sense that the chromosphere and transition region 
are not expected to undergo much damping from this RMHD mechanism; there seems to be insufficient 
time for the turbulent cascade to develop in such small volumes. Dmitruk & Matthaeus (2003) discussed a 
hierarchy of timescales that should be satisfied in order for this turbulent damping process to be strong. A 
simplified version of this hierarchy is as follows; 

to <tR <tf , (52) 



where to is a nonlinear outer-scale eddy cascade time, ta is a timescale for Alfven wave refiection, and tf is 
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representative of the main driving period of the waves. In Figure 14 we compare several of these quantities 
to see where those conditions are met. 

In Figure 14 wc plot a representative reflection time t/j which wc define as 1/|V • V^ij. This is probably 
of the same order of magnitude as the Dmitruk & Matthaeus (2003) "Alfven wave crossing time" tA, since 
the latter can be defined as Ha/Va, or the wave travel time over a representative Alfven-speed scale height 
Ha- Thus, tfl « ^A- For reference we also show the total wave travel time from the photosphere to a given 
height, defined as 



We equate the driving period tf to the spectrum-integrated first-moment period 27r/(a;) (see eq. [41]) and 
we plot its value for the baseline aj = 3 km s^^ model. Interestingly, tf > Ir only below z « 0.5 i?©, and 
the hierarchy is not satisfied at large distances where the damping length is relatively small. Note, though, 
that the demand in the above hierarchy for tf to be the largest of the three timescales comes only from 
the well-known non-WKB result that long periods are the easiest to reflect, and thus they provide the most 
turbulent mixing between inward and outward modes. However, the amount of reflection in our model is a 
known quantity (see, e.g., Figures 7 and 11), and it is irrelevant whether it could have been maximized with 
longer-period waves. Therefore we do not need to consider the last inequality in eq. (52) as a precondition 
for turbulent damping. 

The key piece of the Dmitruk & Matthaeus (2003) timescale hierarchy is that the nonlinear driving time 
to must be short compared to the reflection time t^. Only when this condition is satisfied can the turbulent 
cascade develop to the point where the wave damping can occur efficiently. In Appendix C we describe 
one way of using a model for the anisotropic turbulent power spectrum to derive to (see also Oughton, 
Dmitruk, & Matthaeus 2004). This timescale depends inversely on the outer-scale length and in Figure 
14 we plot the radial dependence of to for two extreme values of the normalizing constant A. An upper 
limit to the outer-scale length is given by A = 1 (i.e., that the "stirring" takes place on the spatial scale 
of the entire network element). A reasonable lower limit for L± at the merging height is the radius of the 
smallest expected MBP flux tube. At the merging height, the flux tubes have expanded to fill the volume 
of the network patch, and thus the MBP radius is given by half the nearest-neighbor distance For the 
lower-limit value of w 350 km (see § 3.1) we obtain a lower limit of A = 0.06. 

We find in Figure 14 that tji begins to exceed to somewhere above a height of 0.1 to 1 Rq. We thus 
solve the nonlinear ly modified wave transport equation only above 0.1 Rq. This restriction allows us to 

make a useful simplifying assumption. As seen in Figure 11, above this height there is only relatively weak 
inward wave power; i.e., {Z-) ^ {Z+). We solve a modified time-steady wave action conservation equation 
in the limit of pure outward propagation. Note also that turbulence does not act in a straightforward 
"monochromatic" way, as was assumed in the constant-frequency linear solutions to eq. (32). We thus ignore 
any frequency (and wavenumber) transport and compute the overall efl:ccts of the turbulent damping on 
the frequency-integrated wave energy (see, however, Verdini & Velli 2003). Thus, we solve the following 
modified spectrum-averaged wave action conservation equation in the limits of no time dependence and 
dominant outward propagation: 




(53) 



uVa d \{u + VAf{Z_) 



(54) 



u + VAdr I uVa 



(see Appendix B). We use the linear solution for {Z+) obtained above in order to solve for the nonlinearly 
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damped value of ^Z-): 



{Z_)l {u + VAfuoVAfi 



exp 




[u{r') + VA{r' 



dr'{Z+{r')) 



)]L^{r') 



(55) 



where all quantities with subscript '0' are given at the effective lower boundary z = 0.1 i?©. This is not a 

completely self-consistent method of obtaining the damped wave power, since no back-reaction on {Z^) is 
computed, but it gives an adequate order-of-magnitude result for the total power in the regions dominated 
by outward propagating waves. 

Figure 15 shows the result of integrating eq. (55) for three choices of A (0.1, 0.35, and 1), as well as 
the solution for no damping (A — > oo). We plot these solutions as velocity amplitudes in order to facilitate 
comparison with the observational data, which are also plotted as in Figure 9. For each value of A, the 
resulting damped C/_ energy density was added to U+ (which remains undamped) to obtain [/tot; and this 
total energy density was split into kinetic and magnetic parts assuming that the linear undamped ratio 
Uk/Ub remains the same for the damped models. The resulting Uk was then used to compute {6V) as in 



The value of A that produces the best agreement with the in situ measurements is approximately 0.35, 
implying a transverse outer scale at the merging height of about 1100 km. This value seems appropriate and 
consistent for motions excited between granules of the same spatial scale. If A were much smaller than this 
value, it would imply that the turbulence was dominated by internal motions in the MBP flux tubes, which 
we essentially ignore in the thin-tube models. Furthermore, note from Figure 15 that for A < 0.1 there is 
so much damping that the approximation ^ {Z+) breaks down above ~100 i?© and we would have 

needed to solve a more complicated set of damping equations. The observational constraints which imply 
A = 0.35 allow us to avoid both of the above consistency problems that arise for smaller values. 

Properties of the A = 0.35 solutions were also plotted in several other figures in order to compare 
with the linear results. Figure 10 shows that the damped magnetic amplitude ratio {SB)/Bq is smaller 
in interplanetary space than the undamped ratio. For the damped model, {SB)/Bq exceeds 1 only above 
135 Rq and does not exceed a value of 1.7 within the computed range of distances (see also Lau & Sireger 
1996). Figure 11 shows the better agreement between the damped Elsasser amplitude {Z-) and the measured 
outward Elsasser energy. The inward amplitude {Z+) still does not agree with the measured inward Elsasser 
energy (see also Hollweg 1990), but note that the linear value was not modified in the damping calculations 
done here. Full simulations of the turbulence would probably see an enhancement of {Z+) and possibly even 
an approach to inward/outward energy equipartition (e.g., Oughton & Matthaeus 1995). 



A self-consistent treatment of waves in the solar atmosphere and solar wind would consider the impact 
of fluctuations on the mean, time-steady plasma properties. Incorporating this back-reaction into the model 
steady state is beyond the scope of this paper, but below we compute a subset of the necessary source terms 
(e.g., the heating rate and the wave pressure acceleration) for use in future models. 



eq. (38). 



6.3. Impact on the Steady-State Plasma 
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6.3.1. Heating due to Turbulent Dissipation 

The nonlinear wave damping discussed in § 6.2.2 produces gradual heat deposition along an open 
magnetic flux tube. The rate of energy conversion is derivable essentially from eq. (54), and here we give a 
slightly more general version of the energy conversion rate based on a comparison between phenomcnological 
turbulence models and numerical simulations. The time rate of change in the general Elsasser variance {Z±)'^ 
is given by 

= -„,(M(2i> (56) 

dt ^ ' 

(e.g., Hossain et al. 1995; Matthaeus et al. 1999; Dmitruk et al. 2001, 2002). Comparison with numerical 
simulations found that the dimensionless parameter a± is of order unity. We thus assume a± = l and define 
the total heating rate Q as 

In our models this quantity is defined consistently only above z = 0.1 -R0, where the damping is applied, but 
below we calculate Q for all heights along the flux tube. We caution that the heating rates below 0.1 Rq 
are only estimates based on the undamped Elsasser amplitudes at those heights. 

Figure 16 plots the heating rate per unit mass (Q/p) for the range of A values shown in Figure 15 in 
addition to an even weaker-damping case of A = 3. For the regions where the nonlinear damping does not 
change the radial proflles of {Z±) from their imdamped values, Q scales with A~^ as expected. However, 
for heights where damping occurs the heating rate is lowered because of the cumulative loss of wave power 
up to that height. Thus, for each height, there is a maximum heating rate that defines a critical value Acrit- 
For A > Acrit the cumulative damping grows weaker, and the decrease in Q comes from the increase in the 
denominator of eq. (57). For A < Acrit, there is more cumulative damping and the decrease in Q comes from 
the loss of local wave power in the numerator of eq. (57). For a distance of r = 2 Rq in the extended corona, 
Acrit is about 0.07, and the maximum heating rate per unit mass at that value is 3.9 x 10^^ erg s~^ g~^. 

We also plot in Figure 16 an empirical heating rate (summed from proton and electron contributions) 

that was derived in order to reproduce a large number of remote-sensing and in situ measurements of fast 
solar wind conditions. We use the heating rates from the SW2 model of Allen, Habbal, & Hu (1998). For 
comparison with the above numbers, the empirically constrained heating rate at r = 2 Rq is about 1.2 x 10^^ 
erg s~^ g~^ (see also Dmitruk et al. 2002; Cranmer 2002). In the extended corona (r « 2-3 Rq), Figure 
16 shows that the empirical heating rate agrees best with the A = 0.35 case that also produces the right 
amount of in situ damping (Figure 15). This rough agreement provides additional support for the outer-scale 
normalization specified by this value. 

It is worthwhile also to compare the heating rates computed above with those computed using a damping 
rate appropriate for isotropic Kolmogorov (1941) turbulence (see, e.g., Hollweg 1986b; Chae et al. 1998; Li 
et al. 1999; Chen & Hu 2001). In order to derive a Kolmogorov heating rate appropriate for quantitative 
comparison with the above results, we make both the quasi-isotropic approximation U- = U+ and we assume 
MHD energy equipartition Uk = Ub- Thus, eq. (57) becomes 

where we use the damped velocity amplitude {5V) as plotted in Figure 15. The above expression differs by a 
factor of \/2 from the simpler formulae given in many other papers; this factor is required in order to remain 



- 28 - 



consistent with the isotropic/equipartition limit of eq. (57). Removing this factor would be equivalent to 
specifying a± = l/-\/2 in eq. (56), but we note that the relative comparison between eqs. (57) and (58) does 
not depend on the value of a±. 

In Figure 17 we plot the comparison between the complete heating rate Q and the Kolmogorov ap- 
proximation Qkoi- For simplicity we show only the case having the preferred A = 0.35 outer-scale normal- 
ization. The curves are substantially different from one another nearly everywhere, which indicates that 
the inward/outward imbalance generated by non-WKB reflection is probably a very important ingredient in 
Alfven wave heating models of the solar wind. The differences are small in the photosphere and low chro- 
mosphere, where strong reflection leads to nearly equal inward and outward wave power. In the extended 
corona, though, the Kolmogorov heating rate begins to exceed the anisotropic turbulent heating rate by as 
much as a factor of 30. The isotropic Kolmogorov assumption assumes the maximal amount of possible mix- 
ing between inward and outward modes, which is inconsistent with the relatively weak reflection computed 
for the corona in our models. 

Figure 17 also shows the radial dependence of an analytic scaling relation for the turbulent heating rate 

given by Dmitruk et al. (2002) for the region above the coronal maximum in Va- This scaling relation gives 
Q oc Bo\dVA/dr\, and we show it for an arbitrary normalization. The analytic heating rate per unit mass 
has a radial gradient less steep than the other curves in Figure 17, but it has a general shape that is similar. 



6.3.2. Wave Pressure Acceleration 

Waves that propagate radially through an inhomogeneous medium exert a nondissipative net force on 
the gas. This net momentum deposition has been studied for several decades for both acoustic and MHD 
waves and is generally called either "wave pressure" or the ponderomotive force (e.g., Bretherton and Garrett 
1968; Dewar 1970; Belcher 1971; Alazraki and Couturier 1971; Jacques 1977). Initial computations of the net 
work done on the bulk fluid have been augmented by calculations of the acceleration imparted to individual 
ion species (e.g., Isenberg & Hollweg 1982; McKenzie 1994; Li et al. 1999), estimates of the departures from 
Maxwellian velocity distributions induced by the waves (Goodrich 1978; Hollweg 1978b), and extensions to 
nonlinearly steepened wave trains (Gail, Ulmschneider, & Cuntz 1990; Pijpers 1995). 

For non-WKB Alfven waves propagating in a radial flux tube, Heinemann & Olbert (1980) gave the 
general second-order expression for the wave pressure acceleration: 

awp = -(v± • V)vi + 7^ [(V X Bi) x B^] (59) 

which has the radial component 

_ 1 d\B^\' f\B^\' d 

In the limit of equipartition between the kinetic and magnetic wave energy densities, the term in parenthe- 
ses is zero. For a spherical geometry (e.g., MacGregor & Charbonneau 1994), we note that the quantity 

i9(ln_Bo)/9r can be written simply as — 2/r. 

In our non-WKB model we have computed the frequency-integrated kinetic and magnetic energy densi- 
ties Uk and Ub. We find that the spectrum- weighted wave pressure acceleration can be expressed in terms 
of these variables as 

p(awp) = -^ + {UB-UK)^{lnBo) . (61) 
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Figure 18 compares the computed values of (flwp) for the undamped {aj = 3 km s~^) and damped (A = 0.35) 
cases, as well as an effective WKB wave pressure computed using just the first term on the right-hand side 
of eq. (61). There is some numerical noise in the plots that comes from the necessity to take numerical 
derivatives. Departures from the idealized WKB form arise only below z w 0.03 Rq and thus may not be 
important for the acceleration of the solar wind. In the photosphere, though, the non-^'^'KB wave pressure 
acceleration exceeds the approximate WKB acceleration by a factor of 100. The WKB approximation also 
gives an unrealistically large spike in (cwp) at the transition region that does not appear in the exact non- 
WKB solution. It remains to be seen if these differences would significantly affect various observational 
constraints placed on the properties of the transition region (e.g., Woods, Holzer, & MacGregor 1990a, 
1990b). 

Figure 18 also shows the magnitude of the Sun's gravitational acceleration ^ as a function of height. 
Presumably, the wave pressure acceleration has little net effect on the corona and solar wind in regions where 
K^wp)! ^ 9- For closed coronal loops, though. Laming (2004) computed a similar non-WKB wave pressure 

acceleration and found that its strength can rival that of gravity even in the chromosphere. The Alfven 
wave reflection in these models is strong enough to produce ion abundance variations that may explain the 
observed dependence on First Ionization Potential (FIP). 



7. SummEiry of Major Results 

The goal of this work was to produce a detailed and self-consistent description of the global energy 
budget of transverse incompressible waves in open magnetic regions of the solar atmosphere. Here we list 
the unique features of the model, key results, and some of the insights gained about the overall wave-plasma 
system. 

1. Measurements of G-band bright points in the photosphere were used to set the power spectrum of 
transverse waves as the lower boundary condition of our model. The observationally inferred power 
spectrum was summed from two phases of MBP motion assumed to be statistically independent: 
isolated random walks and occasional rapid jumps due to MBP merging, fragmenting, or magnetic 

reconncction. 

2. The steady-state plasma density, magnetic field, and flow velocity were constrained empirically from the 
photosphere to a distance of 4 AU. The successive merging of flux tubes on granular and supergranular 
scales in the atmosphere was described using a two-dimensional magnetostatic model of the magnetic 
network element. 

3. Non-WKB wave transport equations, incorporating linear reflection terms, were solved for a grid of 
wave periods between 3 and 3 x 10^ s over the range of heights given above. Below a mid-chromospheric 
merging height the waves were modeled as modified kink-mode fiux tube waves, and above this height 
they were modeled as Alfven waves propagating parallel to the background magnetic field. 

4. The waves are refiected strongly at the transition region, with only about 5% of the wave energy 
transmitted and 95% reflected. Above the transition region, most periods are reflected only weakly by 
the large-scale radial gradient of the Alfven speed, but periods exceeding 1 day are reflected strongly in 
interplanetary space. The period-averaged reflection coefficient from the extended corona (r = 2Rq) 
to 1 AU ranges between 10~* and 10~^. 
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5. At the photosphere, the wave periods containing the most power are between 1 and 30 minutes, but 
the kink-mode waves with periods greater than about 12 minutes are evanescent below the merging 
height. Only the shorter periods propagate to larger distances, and above the merging height the 
power spectrum is dominated by wave periods between about 1 and 6 minutes. In situ observations of 

power spectra that arc strongest at periods of several hours may be explained by many uncorrelated 
flux tubes rotating past spacecraft in the solar wind on these longer timescales. 

6. The period-averaged transverse velocity amplitude of the waves agrees with observed off-limb nonther- 
mal line widths from SUMER and UVCS when the MBP jump amplitude aj is taken to be ~3 km 
s~^. For all reasonable values of aj between and 6 km s~^ the modeled wave amplitudes at heights 
greater than 0.3 AU are significantly larger than the in situ measurements, implying that large-scale 
damping is needed. 

7. We investigated the potential impact of the linear viscous damping of Alfven waves. If the effective 
viscous stress timescale undergoes a transition from its collisionally dominated form to either of the 
possible coUisionless forms (given by the m = 1 or 2 cases of § 6.2.1), the resulting damping lengths 
are always much longer than the local height from the photosphere, and viscosity thus cannot damp 
the waves. 

8. We also considered nonlinear turbulent damping using a phcnomcnological model of the energy loss 
terms in the wave transport equations. The one free parameter is the normalization of the outer-scale 

— 1/2 

correlation length L± of the turbulence, which scales with height as Sq ' . Interestingly, a single choice 
for the constant A = 0.35, which specifies the value of L± at the merging height to be about 1100 km, 

produces both the right amount of damping above 0.3 AU to agree with the in situ measurements and 
the right amount of heating in the extended corona to agree with empirically constrained solar wind 
acceleration models. 

9. Because of the relatively weak degree of reflection in the corona, non-WKB effects are probably not 
important for computing the wave pressure acceleration in regions that contribute to the overall solar 
wind acceleration. However, non-WKB effects produce ordcr-of-magnitude differences (from the WKB 
approximation) in the wave pressure acceleration at and below the transition region. The inclusion of 
damping also affects the total amount of wave pressure acceleration imparted to the solar wind. 

8. Discussion 

The model presented in this paper contains a great deal of the relevant physics of Alfven waves in the 
solar atmosphere and solar wind, but it is incomplete in several ways. In this section we discuss various ways 
to extend the model to improve its physical realism. 

Our cylindrically symmetric magnetostatic model of the network flux tube is highly idealized. A fully 
three-dimensional model of a network "patch," possibly using a high-resolution magnetogram image as the 
lower boundary condition, would have a more realistic geometry (see, e.g., van Ballegooijen & Hasan 2003). 
Such a model would contain a distribution of flux tube strengths, merging heights, and canopy heights. 
The resonant node structure seen in the power spectra (e.g.. Figure 8) would likely disappear in such a 
heterogeneous model. 

Our adopted photospheric boundary condition for transverse flux-tube motions is also an idealized 
approximation. First, the observations used to derive the velocity autocorrelation functions are limited in 
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time resolution. Extremely high-frequency acoustic modes — e.g., with periods below 1 minute — have been 
seen in the photosphere and chromosphere (Deubner 1976; Wunnenberg, Kneer, & Hirzberger 2002; DeForest, 
DePontieu, & Hassler 2003), and the kink- mode spectrum may extend to these frequencies as well. Second, 
much more needs to be learned about the statistics of the "jumps" that we modeled as identical Gaussians. 
A wcll-dcfincd algorithm for measuring the combined power spectrum of both isolated MBP motions and 
their mergings and fragmentings should be developed. It is possible that the transient nature of the motions 
may be better represented by waves having both complex frequency and complex wavenumber (e.g., Wang, 
Ulrich, & Coroniti 1995), or perhaps by suitably defined wavelet functions. Also, new diagnostics such as 
polarization variability within line profiles may add better constraints to our understanding of flux-tube 
motions (Ploner & Solanki 1997). 

Compressible MHD waves and shocks must also be considered alongside the incompressible waves mod- 
eled in this paper (for recent work, see Rosenthal et al. 2002; Bogdan et al. 2002, 2003; Hasan et al. 2003; 
Hasan & Ulmschneider 2004; Wedemeyer et al. 2004; Bloomfield et al. 2004; Suzuki 2004). Longitudinal 
flTix-tubc waves are likely to be excited at or below the photosphere, and these motions can provide addi- 
tional heating (via coUisional damping) and, in some regions, they can transfer some of their energy into new 
Alfvenic fluctuations that were not taken into account in our models. In the solar wind, additional linear 
and nonlinear couplings need to be included to understand better the relatively large measured values of 
{Z+) (see Figure 11). 

The energy and momentum deposition terms computed in § 6.3 need to be included in self-consistent 
solar wind models. Extensive modeling has been done using isotropic Kolmogorov turbulent heating terms 
(effectively, eq. [58]) and the impact of a more realistic treatment of wave reflection is likely to be significant. 
Models including these effects can also be extcinded to the slow solar wind associated with streamers, to 
coronal holes at other phases of the solar cycle (e.g., Miralles, Cranmer, & Kohl 2002), and possibly to 
quasi-steady stnictures like spicules and prominences. 

A substantial unresolved issue is how the wave-related heating and acceleration is apportioned to elec- 
trons, protons, and heavy ions in regions where these species become decoupled. Our original motivation for 
this work was the determination of "outer-scale" (i.e., low- frequency) Alfven wave properties in the corona 
and solar wind, for use as initial conditions in models of MHD turbulent cascade. The kinetic consequences 
of this turbulence have been studied recently by, e.g., Leamon et al. (1999, 2000), Cranmer & van Balle- 
gooijen (2003), Voitenko & Goossens (2003, 2004), Gary & Nishimura (2004), and Gary & Borovsky (2004). 
The ultimate goal of this work is to model the spatial evolution of a reflecting, cascading, and dissipating 
power spectrum of fluctuations as a function of and k±. Observations of preferential ion heating and 
acceleration, as well as departures from isotropic Maxwellian distributions, are key discriminators between 
competing models. 
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of G-band bright point motions that made this work possible. The authors would like to thank Eugene 
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Milano, William Matthaeus, Marco Velli, Wolfgang Kalkofen, S. Peter Gary, and the anonymous referee 
for valuable and inspiring discussions. This work is supported by the National Aeronautics and Space 
Administration (NASA) under grants NAG5-11913, NAG5-12865, NAG5-10996, and NNG04GE84G to the 
Smithsonian Astrophysical Observatory, by Agenzia Spaziale Italiana, and by the Swiss contribution to the 
ESA PRODEX program. This work has made extensive use of NASA's Astrophysics Data System. 
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APPENDICES 




A. Isothermal Kink-Mode Wave Properties 

In this Appendix we discuss a class of analytic solutions to the incompressible thin-tube wave equations 
given in § 5.1. For constant coefficients in eq. (24), we assume that v± depends on time and height as 
giwt-ifcz _ Pqj. nonzero velocity amplitudes, this wave equation becomes a simple quadratic equation in k that 
is equivalent in many ways to the acoustic-gravity wave equation for an isothermal hydrostatic atmosphere 
(e.g., Lamb 1932; Mihalas & Mihalas 1984; Wang et al. 1995). For a real frequency cu, this equation is 
satisfied only for a complex k, with 

k = i^[^±^|-,-l\ (Al) 

where 

"''^^^wt^^V 8^(2^+1) ^^^^ 
is the critical kink-modc cutoff frequency. The second expression above applies for an isolated flux tube in 
an isothermal atmosphere with scale height H = a^/g, where a is the isothermal sound speed and is the 
ratio of gas pressure to magnetic pressure (e.g,. Spruit 1981). In this case, the identity Vph/tOc = yields 
the standard form of eq. (Al). Other useful identities in the isothermal limit, all derivable from transverse 
pressure balance, are 

^ = ^, pVl = 2gH , Bo = V87ra2|Ap| . (A3) 

Only for u) > cJc can there be a real component of k, and thus a propagating wave. A purely imaginary k 
corresponds to evanescence. 

For the backgroimd plasma state described in § 3, Pe/ P = 2.35 and Vph = 6.672 km s"-*^ at the photo- 
sphere. With these values, the critical period (27r/wc) is found to be 12.49 minutes. Above the photosphere, 
the critical period computed with eq. (A2) decreases to a minimum of 9.53 minutes at a height of ~200 km, 
then begins to increase, formally diverging to infinity at the merging height where Ap (see eq. [29]). 
This range of periods compares favorably with the chromospheric estimates of 11.7 minutes given by Spruit 
(1981), derived for /3 = 1, and 9 minutes given by Hasan & Kalkofen (1999), derived for (3 = 0.3. For 
reference to other isothermal models, we use the above relations to derive effective isothermal values of the 
plasma beta and scale height at the photosphere; these are /? = 0.74 and H = 195 km, although the latter 
is about 20% larger than the actual density scale height at z = 0. 

The imaginary part of k gives the time-averaged radial dependence of the velocity amplitude. For 
high frequencies corresponding to propagating kink-mode waves, is proportional to e^^^^ . (Because 
p oc e~^/^, this implies that \v±\ oc as expected from WKB Alfven wave theory.) Thus, the kinetic 

energy density of the waves integrated over the cross-section of the tube (i.e., plv^l"^ / Bq) is constant.'' For 
low frequencies in the evanescent domain, let us define q = {u)/u)c) < 1, and there are two solutions for the 
radial dependence of the amplitude: 

1^;^! cx exp (l ± • (A4) 



^In the isothermal thin-tube Umit, the vertical zero-order magnetic pressure (Bq/Stt) must remain proportional to the local 
gas pressure (pa^), so that Bq oc exp(— z/2i?). 



-33- 



The upper-sign solution is steeper than the propagating solution, and the lower-sign solution is shallower. 
In the limit q ^ 0, the steep solution goes as p"^/^ and the shallow solution is constant. Thus, for the 
steep [shallow] solution, the wave energy density integrated over the tube area grows [decays] with increasing 
height. In § 6 we see that non-WKB low-frequency waves naturally "find" the shallow solution, probably 
because it is the more physical solution with decaying — not diverging — energy density (see also Wang et al. 
1995). 

An additional reason that a stellar atmosphere may "choose" the shallow evanescent solution instead of 
the steep solution was given by Cranmer (1996) for acoustic-gravity waves in a nonmagnetized atmosphere. 
Despite the fact that the solar wind is extremely subsonic in the deep photospheric and chromospheric layers 
we study here, the existence of a nonzero gradient of the outflow speed leads to additional terms in the wave 
dispersion relation. Cranmer (1996) showed that these terms yield a nonzero real part of k for all frequencies, 
even those below the evanescent cutoff. For a subsonic flow {e = u/a <^ 1) in the evanescent regime (g <C 1) 
the real part of the radial wavenumber is given by 

-Seq/H , steep 

+eq/H , shaUow ^ ' 

and the shallow solution is the one with a positive (i.e., upward) phase speed ui/kr corresponding to the 
physically relevant situation of more upward then downward wave power. 

Lastly we compute the frequency-dependent partition between kinetic energy density {Ek — p\v±\'^/2) 
and magnetic energy density {Eb = \B±\'^ /Sir). The linearized induction equation (eq. [27]) is written in 
the isothermal limit as 

u}B_i_ = -kBov± (A6) 

and thus the ratio of energies is 



Ek a;2 

Using the analytic solution for k above, the total energy density is given by 

E„H-i^.Nx{ J^^^^^y^^^, [III (A8) 

where the upper [lower] sign gives the steep [shallow] evanescent solution as defined above. Note that in the 
limit q ^ 0, the shallow solution is dominated by kinetic fluctuations and the steep solution is dominated 
by magnetic fluctuations. 



B. Non-WKB Alfven Wave Transport Equations 

Over the years there have been several versions of the the non-WKB Alfven wave equations published 
in the literature that superficially do not resemble one another. This introduces a potential difliculty in 

comparing the results of these past efforts. In this Appendix we collect several of these differently appearing 
equations and confirm that they arc essentially equivalent. 

The transport equations of Heinemann & Olbert (1980) are expressible in terms of our Elsasser variables 
(eq. [31]) as 

z± = z± [(u ± Ya) ■ V] In (pV4j _ ± v^) • V] In (^p'/^B^'/^) (Bl) 



- + (utV^).V 
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— 1/2 

where we have specified that the transverse length scale of the flux tube is proportional to -Bq . It is 
relatively easy to see the correspondence between the above form and our eq. (32) when we assume that 
the vectors u and Ya are purely radial and that z± is perpendicular to the radial direction. For mass 
conservation in a flux tube, the quantity p^^^B^^^^ is proportional to V^^^^, and thus the above logarithmic 
derivatives yield the scale heights defined in § 5.2. Khabibrakhmanov & Summers (1997) gave essentially 
the same equations, but allowing for more general wave polarization states and pressure anisotropy. 

The Alfvenic transport equations of Zhou & Matthaeus (1990) contain nonlinear terms, and thus in 
some ways are more general than the equations discussed in this paper. However, when looking only at their 
linear terms, one obtains 

v., . v., . (a^) V . (I ± V.) - . (vu± ^) (B.) 

(see also Dmitruk, Milano, & Matthaeus 2001). The above equation is exactly equivalent to eq. (32), but 
some algebraic manipulation is required in order to demonstrate this. Among other steps, one needs to apply 
the auxiliary conditions given by Zhou & Matthaeus (1990), 



u.V(^) = ^V-Q) , (B3) 



together with the fact that V • Bq = 0. 



The non-WKB equations given by VeUi (1993) are similar to those of Zhou & Matthaeus (1990). Re- 
casting Velli's (1993) differently-defined Elsasser variables into our conventional form, one obtains 

dz± . /Zq: — z± 



dt 



+ (utVa)-Vz± = (^^^-^^ V- (|±Va) -z^-V(u±Va) . (B4) 



Only the last term differs from cq. (B2) above, but — at least for the incompressible transverse waves we 
study here — they are equivalent. This can be shown by noting that, for an arbitrary scalar function V(^)) 
one can state that 

z- V(V'B) = tpz-VB (B5) 

only when two conditions are met: (1) the field is divergence-free (V • B = 0), and (2) the fluctuations 
arc purely transverse (i.e., when z • Vtp = 0). In the above comparison, this condition is satisfied because 
ijj = l/y^Anp depends only on the radius r and the incompressible z amplitudes have a zero r-component. 

Barkhudarov (1991) derived the non-WKB transport equations in terms of nonstandard Elsasscr-like 
variables I± = {B± ± v± y/Awp). Also, in that paper the equations were given only for spherical symmetry. 
For a more general superradially expanding flux tube, the equations of Barkhudarov (1991) become 

and these are the speciflc equations solved by the numerical Runge-Kutta code described in § 5.2. 

For completeness, we also give the full equation of wave action conservation (i.e., wave energy conserva- 
tion measured in the comoving wind frame). We use the variable naming convention of Barkhudarov (1991), 
where / and g correspond to outward and inward waves, respectively, but we follow Heinemann & Olbert 
(1980) in defining these quantities in velocity units and for general superradial expansion: 

2 _ \z-\Hu+VAr . _ \z+\'{u-VAr 

^ = uVa ' ^ = uVa • ^^^^ 
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The wave action conservation equation is obtained by writing the four real transport equations (for the real 
and imaginary parts of z+ and Z-) and multiplying each term in the equations by their corresponding real 
or imaginary Elsasser component. The sum of all four equations yields the wave action equation: 



d_ 
di 



pu 



u + Va 



Va 



+ V-[pu{f-g^)\ = 



(B8) 



In regions dominated by outward propagating waves (|z_| » |^+|), where both the background properties 
and wave amplitudes are time-steady, this reduces to simply 

dr 







(B9) 



and the assumption of wave action conservation is equivalent to the assumption that / is constant (see also 
Jacques 1977; Barkhudarov 1991). This is not necessarily a WKB result, but it applies only in regions of 
negligible reflection. 



C. Issues Relating to Anisotropic MHD Turbulence 

In the wave transport models presented in this paper we have not been overly concerned with the 
direction of the Alfven wavevector k relative to the background magnetic field Bq. Mainly this is because 
the dispersion relation of the waves in the MHD limit depends only on , the parallel component of the 
wavevector, 

u = {uTVA)k\\ (CI) 

for the z± modes, and not on the component k± perpendicular to the field. Thus the monochromatic wave 
transport equations "know" about fcy via their dependence on wave frequency u> but are independent of k± . 
In this Appendix we investigate several issues concerning the expected distribution of propagation directions 
and its evolution as MHD turbulent cascade develops. 

Prom a linear standpoint, an initial distribution of propagation directions 6 (i.e., the angle between k 

and Bq) will evolve with heliocentric distance due to several physical processes. The radial variations of 
the inertial-frame phase speeds u ^ Va induce a radial "stretching" in wavelength and thus a large-scale 
decrease in fc|| with increasing height. The transverse spreading of the open flux tube should also cause the 
perpendicular wavenumber k± to decrease with increasing height as well. The overall sense of increase or 
decrease in depends on the time-steady properties of the magnetic field and outfiow. Other effects can 
contribute to a radial variation in 9. First, any transverse gradients of the Alfven speed can refract oblique 
waves toward the center of the flux tube, thus decreasing 9 (Wentzel 1989). Second, viscous damping of 
high-frequency waves preferentially dissipates oblique fluctuations, since L^. oc cos^ 9. The analysis in § 6.2.1 
assumed 9 = 0, but the impact of viscous damping can be enhanced significantly if the waves arc sufficiently 
oblique. All of the above effects should be included if the radial variation of 6 (or (9) averaged over the 
power spectrum) is to be computed. 

As the turbulent cascade develops in the corona, an additional nonlinear evolution of the wave energy 
within "wavenumber space" {k\\,k±) occurs. It has been known for several decades that MHD turbulence in 
the presence of a steady-state magnetic field develops a strong anisotropy in wavenumber. Both numerical 
simulations and RMHD analytic descriptions indicate that the spectral transport from large to small spatial 
scales proceeds mainly in the two-dimensional plane perpendicular to Bq (see, e.g., Higdon 1984; Shebalin, 
Matthaeus, & Montgomery 1983; Goldreich & Sridhar 1995, 1997; Bhattacharjee & Ng 2001; Cho, Lazarian, 
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& Vishniac 2002; Oughton et al. 2004). That is, the cascade spreads out the power to successively larger 
values of kj_ while leaving /cy relatively unchanged. 

Here wc follow Cranmcr & van Ballcgooijcn (2003) and give an analytic solution for the turbulent power 
spectrum in the limit that the cascade is allowed to proceed to its final "driven" steady state in a small 
homogeneous volume of plasma. We then compare the frequency dependence of this power spectrum to 
the empirically constrained (linear) frequency spectrum as presented in Figure 8. Let us define the three- 
dimensional total power spectrum as 

Utot = J c;'kP3D(k) (C2) 

where we write the volume clement d'^k in cylindrical coordinates as 2irk± dk±dk\\ and thus assume symmetry 
in the two directions transverse to the background field. Cranmer & Ballegooijen (2003), following the general 
development of RMHD anisotropic cascade theory, assumed that the turbulence eventually leads to a state 
of complete mixing between the inward and outward modes, and thus (7+ = U- = t/tot/2. This assumption 
is incompatible with the non-WKB wave reflection models computed in this paper, but wc follow the earlier 
analysis in order to compare the shapes of the relevant power spectra and to obtain insight about how to 
better model the turbulent cascade. 

For fully developed anisotropic MHD turbulence, we assume the power spectrum to be a separable 
function of two variables: k± and a nonlinearity parameter y defined as the ratio of the local wind-frame 
frequency V^fcii to an assumed nonlinear eddy turnover rate {SV)k± (see, e.g., Goldreich & Sridhar 1995). 
We use the notation of § 2.3 of Cranmer & van Ballegooijen (2003) and define 

P3D(fc||,fc^) = ^^^^^g{y) (C3) 
where W±{k±) is a reduced power spectrum, scaled to velocity-squared units and defined as 

1.2 r + oo 

W±(fc±) = / rf/C||P3D(fc||,fci) . (C4) 

For the MHD inertial range, W±_ is proportional to kj^^^. For simplicity, let us define 

W^{k^) = I ll^^^kJko.,)-y^{3np)-^ , fci > kZ ^^^^ 

for a given outer-scale perpendicular wavenumber fcout- The factor of Swp above is needed to normalize the 
full power spectrum as defined in eq. (C2). 

The fc|| dependence of the power spectrum is contained in the dimensionlcss g{y) function in eq. (C3). 
The condition y = 1 is defined as "critical balance" by Goldreich & Sridhar (1995), and their analysis only 
constrains the general form of g{y), not its exact value. The bulk of the wave power is believed to reside at 
low values of y (e.g., g is nonzero only for \y\ < 1), where we define 

This condition captures the highly nonlinear state of turbulence, for which a coherent wave survives for no 
more than about one period before nonlinear processes transfer its energy to smaller scales. One reasonable 
possibility for the parametric dependence of this function is a normalized Gaussian, 

g{y) « ^ e-y' . (C7) 
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Cho et al. (2002) analyzed numerical simulations of MHD turbulence and found that g{y) can be fit reasonably 
well by decaying exponential or Castaing functions. Cranmer & van Ballegooijen (2003) solved a simple 
wavenumber diffusion equation to obtain the analytic form 



9{y) = 



r(n) 



r(n- 0.5) 



7(1 - cy 



1 + 



7(1 -C)V 



a 



which is normalized to unity when integrated over all y, and with 

_ (/3/7) + C + 1 
2(1-0 ■ 



(C8) 



(C9) 



For the MHD inertial range, ( = 1/3, and the three constants a, (3, and 7 describe the relative strengths of 
fc|| diffusion, k± advection, and k± diffusion, respectively. We assume a = 7, thus leaving the only "free" 
parameter to be the ratio P/j. Cranmer & van Ballegooijen (2003) discussed the most realistic values of 
this ratio; earlier models implied that f3/j « 1, but one would need this ratio to be smaller than about 
0.25 in order to produce enough parallel cascade in the corona to heat protons and heavy ions via cyclotron 
resonance. 

With the above definitions it becomes possible to derive the effective frequency power spectra P±{u}) 
that are consistent with fully developed anisotropic turbulence. Noting that we define w as the frequency 
measured in the inertial frame centered on the Sun, we define 



P±{<^) = \j rf'kP3D(k) 5 [w - {uTVA)k\{\ 



(CIO) 



where the Dirac delta function collapses the three-dimensional wavenumber spectrum into a one-dimensional 
frequency spectrum. This "translation" must be done differently for the inward and outward mode spectra, as 
indicated by the difference between the inertial-frame phase speeds of the two modes. The above expression 
is defined formally for both positive and negative frequencies; below we multiply by 2 in order to consider 
only positive frequencies. For the purposes of this Appendix, let us use the simpler form of eq. (C7) for g{y) 
that gives the straightforward analytic solution 



P±{oj) = 



exp 



'i'out.it 



where we define the effective outer-scale frequency 



Va 




(Cll) 



(C12) 



and wc assume fcout = I-k/L^. We note that the use of cq. (C8) for g{y) would have yielded a solution similar 
to the above form of P±, only differing slightly in a small region of frequency near Wout.i- This frequency 
sets the scale for the spectrum. We use the effective outer-scale frequency for outward propagating waves to 
define the nonlinear "driving time" to = 2-jT/u)out,- plotted in Figure 14. 

Figure 19 compares the above turbulent frequency power spectra P±{u)) with the empirically constrained 
linear frequency spectra P± (cj) as derived in § 5 above. For simplicity we plot the spectra at only one 
heliocentric distance (r = 2Rq) and note that one can use the radial dependences of the various timescales 
plotted in Figure 14 in order to extend these results to other heights. The above assumption of full mixing 
for the nonlinear cascade (i.e., U- = U^) is far from the linear non-WKB result of U- f/+ in the extended 
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corona. The inability of the analytic RMHD anisotropic cascade theory to take account of the strong 
imbalance between inward and outward wave amplitudes is a major shortcoming of this kind of model. We 
need a better understanding of the detailed wavenumber dependence of the inward-to-outward power ratio 
for MHD turbulence. 

In Figure 19 we use the preferred normalization A = 0.35 to set the scale of L±. With this normalization 
at r = 2Rq, the outer-scale driving timescales 27r/wout,± are 5.7 and 4.1 minutes for the + and — signs, 
respectively. The peaks of the scaled power spectra ujP± occur at frequencies of ^ 1.12a;out,±7 but the 
empirically constrained linear spectra loP± contain most of their power at slightly higher frequencies (i.e., 
shorter periods). The fact that the two spectra exhibit maxima within the same order of magnitude of 
frequency, though, may be just a result of choosing the radius of r = 2Rq. Figure 14 shows that if the 
curves in Figure 19 were computed instead at 1 AU, the peaks of the linear P± curves would be in the same 
range of periods (1 to 5 minutes), but the peaks of the P± curves would be at periods between 30 and 40 
minutes. At 1 AU the bulk of the power in the linear frequency spectrum would be at an extremely high 
local frequency from the standpoint of the cascade (i.e., y ^ 1). This suggests that the Sun emits a "fossil" 
frequency spectrum that eventually can advect into a region of the three-dimensional wavenumber spectrum 
that was previously believed to be strongly depleted due to the anisotropic cascade. The implications of this 
result will be investigated in future work. 
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i ~1.5 Mm y i ~30 Mm Y i ~5000 Mm )- 

Fig. 1. — Summary sketch of the magnetic field structure discussed in this paper, with the fields of view successively 
widening from (a) to (b) to (c). (a) Intergranular lanes host G-band bright points that are shaken transversely 
to generate kink-mode waves, (b) Above the height where individual flux tubes merge, the coronal-hole network 
field is mainly open (with a funnel/canopy structure), and kink- mode waves are transformed into Alfvon waves, (c) 
Non-WKB waves in the solar wind propagate and reflect depending on their frequencies, and MHD turbulent cascade 
can occur where outward and inward waves interact nonlinearly. The inverted solar image was obtained by EIT 
(Extreme-ultraviolet Imaging Telescope) on SOHO (e.g., Moses et al. 1997). Field lines in (c) are plotted at 2.1° 
intervals at the solar surface, and thus each pair of lines encompasses -^1-2 network funnels. 
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Fig. 2. — Magnetic field lines in a cross section of the magnetostatic model of a network element. Dotted lines at 
z = 0.6 and 1 Mm denote the merging height Zm and the canopy height Zc- Above Zm the model has cylindrical 
symmetry about the left axis of the plot; below Zm we model the flux tubes as evenly distributed throughout the 
circular network patch of radius ~3 Mm, and the plotted cross sections are shown only for illustrative purposes. 
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Fig. 3. — Steady-state plasma conditions along the central axis of the network flux tube, (a) Magnetic field strength, 
multiplied by x^, for the adopted model {solid line), the model of Hackenberg et al. (2000) {dotted line), and the model 
of Li (2003) {dashed line), (b) Adopted model of the hydrogen number density {solid line), empirical limits taken 
from the minima and maxima of observations by Guhathakurta & Holzer (1994), Fisher & Guhathakurta (1995), 
and Doyle et al. (1999) {dotted lines), and the quantity ptot/l-2/mH {dashed line), (c) Alfven speed Va {solid line), 
outflow speed u {dotted line), and flux tube phase speed Vph {dashed line). 
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Fig. 4. — Outline of the empirical procedure used to derive the photospheric power spectrum of MBP kink-mode 
wave energy. Horizontal MBP positions (a) are differentiated to obtain velocities (b), which axe separated into "walk" 
and "jump" phases (c,d). The autocorrelation functions (e,f) of these velocity time series are computed and Fourier 
transformed to obtain power spectra (g,h) as a function of frequency lu, here plotted versus period {2n /uj) in minutes. 
The kinetic energy power spectra from the statistically independent walk {dashed lines) and jump {dotted lines) 
phases are summed (i) and used to compute the total energy spectrum {solid line). 
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Fig. 5. - Height dependence of the outward propagating wave action flux for a selection of periods, each normalized 
to its value in the photosphere. (See labels for values of the wave period.) The merging height Zm and the transition 
region (T.R.) are labeled with arrows. 




Fig. 6. — Contour/grayscale plot of the ratio Ek/Eb as a function of height and wave period. The MHD-like regime 

of energy equipartition (i.e., a ratio of nearly 1) is shown in light gray, bounded by solid contours at values of 0.8 and 
1.2. The kinetic-dominated regime is in white with dotted contours at values of 10^ and IC*. The magnetic-dominated 
regime is in dark gray with dashed contours at a value of 0.1. 
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Fig. 7. — Contour/grayscale plot of the ratio E+/E- as a function of heigiit and wave period. The region of "strong" 
reflection is in gray, bounded by solid contours at values of 0.5 and 0.75. The location of the maximum value (0.99825) 
is denoted by a star. Other contours range between values of 10~* and 10~^ and are labeled by their values. The 
transition region is clearly seen at « « 0.003 i?©. 




Fig. 8. — Period and height dependence of the scaled total power spectrum u]Ptot{(^,r) for six heights (top to 
bottom): z = 0, 0.6 Mm (merging height), 2.15 Mm (transition region), 1 Rq, 9 Rq, and 214 Rq. The latter 3 
heights are labeled with their heliocentric radii. The MBP jump amplitude was assumed to be = 3 km s~^. 
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Fig. 9. — Height dependence of the frequency-integrated velocity amplitude. Solid lines give the undamped value 
of (SV) for three choices of aj-. 0, 3, and 6 km (from bottom to top). The dashed line is ((5V)s for the aj = 3 
km s~^ case. Other lines and symbols correspond to observations discussed in detail in § 6.1. 




Fig. 10. — Frequency-integrated ratios for the baseline aj = 3 km s^^ model: kinetic-to-magnetic energy density ratio 
{dotted line), dimensionless magnetic amplitude for the undamped linear model {solid line) and for the nonlinearly 
damped A = 0.35 model {dashed line). The merging height, transition region and 1 AU are labeled with thin vertical 
lines. 
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Fig. 11. — Height dependence of the frequency- integrated Elsasser variables for the CTj = 3 km s"'^ model: linear 

undamped {Z-) {solid line), nonlinoarly damped {Z-) {dashed line), WKB height dependence of Z- {dot-dashed 
line), and the linear value of {Z+) {dotted line). In situ measurements from Helios and Ulysses are shown as solid 
bars (Bavassano et al. 2000), with the upper set corresponding to Z- and the lower set to Z+. 
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Fig. 12. — Energy flux density of frequency-integrated waves for the = 3 km s~^ model: flux within the MBP tubes 

{solid line), averaged over granular spatial scales {dashed line), and averaged over the supergranule funnel/canopy 
structure {dotted line). The slope discontinuities at Zm (i.e., v/Rq — 1 = 8.6 x 10"'*) are due to the discontinuity in 
dBo/dr at the merging height. 
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Fig. 13. — Viscous damping length scales plotted versus height for three assumptions about the effective viscous 
timescale: classical m = viscosity {dot-dashed line), Williams' (1995) m = 1 viscosity {dashed line), and our 
transverse diffusion m = 2 viscosity {solid line). Also plotted is the local height z {dotted line). 
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Fig. 14. — Comparison of approximate timescales important for nonlinear damping: tn is the wave reflection time 
{thick solid line), tf is the spectrum-weighted wave period {thin solid line), ttrav is the total wave travel time from the 
photosphere {dashed line), and to {gray region) is a nonlinear turbulent driving time computed for A = 0.06 {bottom 
dotted line) and A = 1 {top dotted line). 
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Fig. 15. — Height dependence of the frequency-integrated velocity ampUtude for aj = 3 km and a range of values 
of the dimensionless outer-scale length-scale normalization constant A: 0.1 (dashed line), 0.35 [dash-dotted line), 1 
{dash-triple- dotted line), and no damping {solid line). Other lines and symbols correspond to observations discussed 
in detail in § 6.1. 




Fig. 16. — Turbulent heating rate per unit mass for the models shown in Figure 15. The same line-styles are used 
for A = 0.1, 0.35, and 1, and we also show the weaker damping case A = 3 {solid line). The diamonds illustrate the 
empirically constrained heating rate used by Allen et al. (1998). 
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Fig. 17. — Turbulent heating rate per unit mass for the A = 0.35 model {solid line) and for the analogous isotropie 
Kolmogorov heating rate {dashed line). The radial dependence of the Dmitruk et al. (2002) analytic formula is also 
shown with an arbitrary normalization {dotted line). 




Fig. 18. — Wave pressure acceleration versus height for the undamped aj — 3 km s~^ model {solid line), for the 
A = 0.35 damped model {dashed line), and for the undamped model computed using the WKB approximation for 
the acceleration {dot-dashed line). The magnitude of the Sun's gravitational acceleration is also plotted {dotted line). 
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Fig. 19. Comparison of scaled fluctuation power spectra at 2 Rq. The linear, empirically constrained frequency 
spectra uiP- {solid line) and ujP+ (dashed line) are plotted alongside the nonlinear cascade frequency spectra projected 
from 3D wavenumber space at the same height: ujP- (dotted line) and ujP+ (dash-dotted line). 



